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1.  INTRODUCTION 


Fatigue  is  a  process  by  which  a  component  is  weakened  by  repeated 
loading.  The  stress  can  be  much  below  the  ultimate  strength  or  even  the 
yield  stress  of  the  material.  However,  over  a  large  number  of  cycles  the 
results  can  be  catastrophic.  Fatigue  has  to  be  considered  as  an  important 
problem  based  on  the  economy  and  safety  of  components  and  structures. 
Effective  control  of  the  fatigue  failures  starts  from  the  design  phase,  the 
material  and  the  processes  to  be  used  for  manufacturing  a  component 
should  aim  towards  the  fatigue  tolerant  design.  The  fatigue  failure  still 
might  occur  in  a  well  designed  component,  reasons  for  which  can  be 
severe  utilization,  accidental  loads,  aggressive  environment,  poor 
maintenance  and  inspection.  As  a  result  the  user’s  fatigue  problems  are 
more  associated  with  maintenance,  inspections  and  non-destructive  testing 
techniques.  The  main  concerns  for  the  user  can  be:  where  the  crack  will 

originate,  how  fast  will  the  crack  grow  and  what  is  the  critical  size  of  the 
crack? 

Forest  [1]  in  his  studies  showed  that  crack  nucleation  starts  early  in  the 
fatigue  life  but  most  of  the  fatigue  life  is  spent  in  the  micro  range.  For 
convenience  sake  it  is  termed  as  the  nucleation  phase  or  period  of  the 
micro-crack  growth. 

The  different  phases  of  growth  of  fatigue  crack  from  inception  to  final 
fracture  are  illustrated  in  the  Fig.l. 
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Figure  1.  Different  phases  of  the  fatigue  crack  growth  [2]. 

The  crack  is  termed  as  a  macro-crack  when  it  is  large  enough  to  be  seen 
by  the  naked  eye  according  to  the  classical  definition.  According  to 
another  definition  the  crack  is  considered  as  a  macro-crack  when  it  has 
sufficient  depth  or  length  at  which  the  local  conditions  for  the  crack 
nucleation  no  longer  affect  the  crack  growth.  In  other  words  the  crack 
growth  will  depend  on  the  bulk  properties  of  the  material.  Another 
definition  which  is  very  important  to  the  context  of  the  study  is  that  “a 
crack  is  a  macro-crack  as  soon  as  fracture  mechanics  is  applicable”.  It  is 
rephrased  as  “A  crack  is  a  macro-crack  as  soon  as  the  stress  intensity 
factor  K  has  a  real  meaning  for  describing  its  growth”  Schijve  [2], 
Depending  upon  the  design  and  the  function  of  the  component,  the 
percentage  of  the  fatigue  life  spent  in  different  phases  may  vary.  A  simple 
representation  of  few  such  cases  is  shown  in  the  Fig. 2. 
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Figure  2.  A  simple  representation  of  different  crack  growth  cases 


1.1  CHARACTERISTICS  OF  FATIGUE 

The  following  characteristics  are  common  to  fatigue  in  all  materials: 

(a)  The  fatigue  process  starts  with  a  microscopic  crack,  called  the 

initiation  site,  which  then  grows  with  each  subsequent  movement, 

,) 

which  can  be  analyzed  by  fracture  mechanics. 

(b)  Fatigue  crack  grows  as  a  result  of  cumulative  damage,  and  the 
material  does  not  recover  when  rested.  The  greater  the  applied 
stress,  the  shorter  the  life. 

(c)  Fatigue  life  is  affected  by  various  factors  such  as  temperature, 
environment,  mean  stress,  loading  mode,  surface  finish  etc. 


(d)  Some  materials  like  steel  and  titanium  exhibit  a  fatigue  limit,  a 

stress  below  which  there  is  no  fatigue  crack  growth  due  to  repeated 
loading. 

1.2  PHYSICAL  NATURE  OF  FATIGUE  DAMAGE 

At  microstructure  level  all  materials  are  anisotropic.  Within  the  grains 
there  is  anisotropy  along  the  crystal  planes.  So  when  viewed  at 
sufficiently  small  size  the  non  uniform  microstructure  or  inclusions  result 
in  non-uniform  distribution  of  stresses.  Regions  where  the  stresses  are 
severe  are  usually  the  points  where  the  fatigue  damage  starts. 

For  ductile  engineering  metals,  grains  that  have  unfavorable  orientation 
relative  to  the  applied  stress  first  develop  the  slip  bands.  Individual  slip 
bands  become  more  severe  and  some  develop  into  cracks  within  grains, 
which  spread  to  other  grains,  joining  with  other  similar  cracks  to  grow 
into  a  large  crack  which  can  propagate  to  failure. 

For  the  hard  and  brittle  materials  with  somewhat  limited  ductility,  the 
damage  due  to  slip  bands  or  within  grains  is  very  less.  Damage  is  more 
concentrated  around  the  defects  in  the  material.  A  small  crack  can  develop 
at  a  void,  inclusion,  grain  boundary  or  a  surface  scratch  and  generally 
grows  in  a  plane  normal  to  the  applied  stress. 
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1.3  FEATURES  OF  A  FATIGUE  CRACKED  SURFACE 


Where  a  failure  is  dominated  by  growth  of  crack,  the  resulting  fracture 
surface  when  viewed  microscopically,  generally  exhibits  a  smooth  area 
near  its  origin.  The  crack  growth  region  is  fairly  flat  and  is  oriented 
normal  to  the  applied  tensile  stress.  Rougher  surfaces  generally  indicate  a 
more  rapid  growth. 


c 

O) 
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Figure  3.  Fracture  surface  showing  fatigue  crack  growth  and  final 
brittle  failure  in  18  Mn.  Steel  [3]. 

Curved  lines  concentric  about  the  crack  origin,  called  ‘beach  marks’  are 
often  present  and  mark  the  progress  of  crack  at  various  stages.  Beach 
marks  indicate  changes  in  the  texture  of  fracture  surface  as  a  result  of 
crack  being  accelerated  or  retarded,  which  may  occur  due  to  an  altered 
stress  level,  altered  temperature,  environment  or  due  to  some  other  cause. 


5 


Beach  marks 


Figure  4(a).  Beach  marks  and  the  crack  initiation  sites  on  a 
fractured  surface  [4] 

If  the  fracture  surface  in  a  ductile  material  is  viewed  at  a  microscopic 
level,  it  often  reveals  the  presence  of  marks  left  by  the  progress  of  crack 
with  each  cycle.  These  marks  are  called  ‘striations’  and  can  be  seen  in  the 
Fig.4(b). 


Striations 
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Figure  4(b).  Fatigue  striations  on  a  highly  magnified  fracture  surface  [4] 
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1.4  STRESS  FIELDS  AHEAD  OF  CRACK  TIP 

If  a  crack  is  considered  to  exist  in  a  material,  and  a  uniform  stress,  S,  is 
applied  remotely.  The  uniform  stress  field  is  altered  in  the  vicinity  of  the 
crack;  stress  rises  sharply  near  the  crack  and  has  a  maximum  value  close 
to  the  crack  tip.  A  sharp  crack  causes  a  severe  concentration  of  stress  and 
is  theoretically  infinite  if  the  crack  is  ideally  sharp. 

From  the  theory  of  fracture  mechanics,  a  quantity  called  “stress  intensity 
factor”,  K,  can  be  defined  that  characterizes  the  severity  of  crack 
situation  as  affected  by  crack  size,  stress  and  geometry.  In  defining  K,  the 
material  is  assumed  to  behave  in  a  linear-elastic  manner  according  to 
Hooke’s  law,  and  the  approach  used  is  called  “linear  elastic  fracture 
mechanics”,  LEFM.  K  is  given  by: 

K  =  FSjlTia.  (la) 

where  F  is  a  geometry  factor,  F=1  for  an  infinite  plate  and  ‘a’  is  the  crack 
length. 

However,  an  infinite  stress  cannot  exist  in  a  real  material.  In  ductile 
materials  plastic  deformations  occur  in  the  vicinity  of  the  crack  tip.  The 
region  within  which  the  material  yields  is  called  the  plastic  zone.  Due  to 
this  yielding  and  blunting  of  the  crack  tip,  the  infinite  stress  is  reduced  to 
a  finite  value  close  to  the  crack  tip. 

In  a  loading  part  of  the  cycle  the  distance  ahead  of  the  crack  where 
material  yields  is  called  the  monotonic  plastic  zone,  rm.  During  the 
unloading  part  of  the  load  cycle  the  elastic  material  around  the  monotonic 
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plastic  zone  pushes  on  to  the  already  yielded  material,  and  a  small 
reversed  cyclic  zone  is  formed  ahead  of  the  crack  due  to  material  yielding 
in  compression.  This  zone  is  called  the  cyclic  or  reversed  plastic  zone,  rc. 
An  illustration  of  the  typical  elastic  stress  fields  ahead  of  the  crack  and 
the  formation  of  the  plastic  zone  are  shown  in  the  Fig. 5(a)  ahead. 


Figure  5(a).  Illustration  of  the  elastic  stress  field  ahead  of  the  crack, 
and  the  redistribution  of  stress  due  to  plastic  zone  formation. 

The  formation  of  monotonic  and  cyclic  plastic  zones  ahead  of  the  crack 
tip  is  shown  in  Fig. 5(b).  Also,  a  small  region  exists  ahead  of  the  crack  tip 
where  due  to  crack  tip  blunting  the  stress  and  strain  are  approximately 
constant.  This  small  region  is  called  process  zone,  p*. 
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Monotonic  plastic  zone,  u 


Figure  5(b).  Different  zones  ahead  of  the  crack  tip  under  cyclic  load. 

1.5  MEASUREMENT  OF  FATIGUE  BEHAVIOR  OF  MATERIALS 
To  measure  the  fatigue  resistance  of  materials,  a  cyclic  load  is  applied  to 
a  specimen  until  it  breaks.  The  number  of  cycles  to  failure,  called  the 
fatigue  life,  Nf,  is  recorded.  The  logarithm  of  this  life  is  plotted  against 
the  stress  (or  the  log  of  stress)  to  develop  an  S-N  curve  (stress  versus 
number  of  cycles  to  failure)  as  shown  in  Fig. 6.  This  curve  defines  the 
fatigue  resistance  of  a  material. 
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Figure  6.  S-N  curve  for  an  unnotched  A517  steel  specimen  [5] 

As  can  be  seen  in  Fig. 6,  there  is  a  great  deal  of  scatter  in  the  fatigue  life. 
Consequently,  a  large  number  of  tests  have  to  be  run  in  order  to  determine 
the  fatigue  resistance. 

For  some  materials,  such  as  steel,  there  is  a  stress  below  which  no  fatigue 
failure  is  observed.  This  is  called  the  ‘endurance  stress’  or  the  ‘endurance 
limit’,  in  terms  of  the  stress  or  the  number  of  cycles  to  failure, 
respectively.  The  endurance  limit  is  typically  on  the  order  of  10^  cycles. 
Other  nonferrous  materials,  such  as  A1  and  Al-alloys,  do  not  have  a  true 
endurance  limit.  In  other  words,  the  S-N  curve  is  never  horizontal.  For 
these  materials,  an  artificial  endurance  limit  is  defined  as  the  stress  that 
will  cause  fatigue  failure  in  a  given  number  of  cycles,  usuallylO’  or  10* 
cycles. 
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1.6  CRACK  GROWTH  RATE  CURVE 


Fracture  mechanics  is  important  in  life  prediction  of  components  that  are 
subjected  to  time  dependent  crack  growth  mechanism  such  as  fatigue  or 
stress  corrosion  cracking. 

Crack  growth  can  be  caused  by  a  cyclic  loading  behavior  called  fatigue 
crack  growth.  The  rate  of  cracking  is  an  important  parameter  to  quantify 
the  time  required  before  the  material  will  fail.  The  rate  of  cracking  is 
usually  called  the  crack  growth  rate,  da/dN,  represents  the  crack  extension 
per  cycle. 

Paris  et  al.  [6]  proposed  a  relation  to  describe  crack  growth  in  metals  with 
respect  to  the  applied  stress  intensity  range  of  the  following  type: 

da/dN  =  C(AK)™  (lb) 

where  da/dN  is  the  crack  growth  per  cycle,  AK  is  the  SIF  range;  C  and  m 
are  the  material  constants. 

The  crack  growth  behavior  for  a  material  can  be  described  on  a  log-log 
plot  as  shown  in  the  Fig. 7. 
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da/dN 


Figure  7.  Illustration  of  fatigue  crack  growth  rate  curve. 

At  low  growth  rates  i.e.  in  the  near  threshold  region  the  curve  generally 
becomes  steep  and  appears  to  approach  a  vertical  asymptote,  denoted  by 
AKth  which  is  called  the  threshold  SIF  range.  AKth  is  usually  interpreted 
as  a  lower  limiting  value  of  AK  below  which  the  crack  growth  is  assumed 
to  be  zero. 

In  the  intermediate  region  of  the  curve,  the  crack  growth  rate  follows 
almost  a  straight  line  plot  in  log-log  scale.  This  region  is  called  the  ‘Paris 
Region’  where  there  is  a  stable  crack  growth  and  can  be  defined  by  Eq.l. 
At  high  crack  growth  rates  above  the  Paris  region  there  is  rapid  unstable 
cracking  and  the  curve  again  becomes  steep  just  prior  to  the  final  failure 
when  the  material  reaches  the  critical  value  of  K  called  the  “fracture 
toughness”,  Kc.  This  is  called  the  static  failure  region  of  the  curve. 
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2.  LITERATURE  REVIEW 


In  order  to  have  damage  tolerant  design,  it  is  important  to  know  the  crack 
growth  behavior  of  a  material,  so  that  the  component  can  be  changed  or 
repaired  before  the  crack  reaches  a  critical  size.  For  this,  it  is  necessary 
to  determine  the  rate  of  crack  growth  with  as  much  reliability  and 
accuracy  as  possible.  During  the  last  thirty  years  or  so,  considerable 
attention  has  been  devoted  to  investigations  into  crack  propagation 
relationships  under  cyclic  loading.  Most  of  these  studies  attempt  to  relate 
the  crack  propagation  rates  with  some  function  of  stress  intensity  factor 
range,  AK,  based  upon  experimental  data.  Although  the  models  are 
applicable  for  a  particular  set  of  conditions  for  which  they  were  derived, 
they  lack  generality  and  require  experimental  determination  of 
coefficients  and  exponents  which  are  not  related  to  the  basic  material 
properties.  There  are  many  crack  growth  softwares  available,  which  use 
different  fatigue  crack  growth  models  to  predict  the  failure  life,  but  none 
has  the  desired  predicting  accuracy.  Either  these  require  extensive 
experimental  data  or  some  kind  of  special  fitting,  coupled  with  lots  of 
fatigue  experience  to  predict  the  crack  growth  behavior.  To  develop  a 
better  predictive  model  and  to  shed  some  light  on  the  damage 
accumulation  process,  various  fatigue  crack  propagation  models  have  been 
studied,  a  review  of  these  models  will  be  presented  in  the  following 
discussion. 
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Fatigue  crack  growth  is  a  phenomenon  of  damage  and  separation  of  the 
material  ahead  of  a  crack  under  cyclic  loading.  It  has  been  reported  [7,8] 
that  plastic  strains  are  the  cause  of  damage  accumulation  at  all  levels  of 
loading,  and  that  the  significant  damage  takes  place  when  the  plastic 
strain  exceeds  a  threshold  value  for  the  material  enveloped  by  the 
reversed  or  cyclic  plastic  zone.  Many  researchers  [7-20]  support  the  fact 
that  the  greatest  damage  in  the  material  is  done  in  a  small  region  ahead  of 
the  crack  tip,  called  the  ‘process  zone’.  The  material  in  this  zone 
undergoes  high  plastic  deformation  that  damages  the  material,  As  soon  as 
this  accumulated  damage  reaches  a  critical  level,  the  material  element 
fails. 

In  the  case  of  small  scale  yielding  i.e.  where  the  plastic  zone  ahead  of  the 
crack  tip  is  small  as  compared  to  the  component  size,  represents  most 
crack  growth  situations.  Plastic  zone  is  surrounded  by  a  much  larger 
volume  of  elastically  deformed  material.  The  plastic  strains  in  the 
monotonic  and  reversed  plastic  zones  are  controlled  by  the  surrounding 
elastic  deformation  fields,  which  in  turn  are  related  to  the  stress  intensity 
factor,  K. 

Rice  [7]  traced  the  deformation  history  of  a  material  element  after  it 
enters  the  reversed  plastic  zone.  Dugdale  model  was  used  to  calculate  the 
total  absorbed  plastic  strain  energy,  and  it  was  assumed  that  the  fatigue 
damage  is  proportional  to  this  energy.  He  obtained  a  fatigue  crack  growth 
rate,  da/dN,  proportional  to  fourth  power  of  AK.  Due  to  lack  of  fully 
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elastic-plastic  solutions  of  crack  tip  fields  for  Mode  I  loading  various 
assumptions  of  strain  distributions  were  made  in  this  model.  It  was  also 
assumed  that  the  deformation  fields  ahead  of  the  stationary  crack  and  a 
propagating  fatigue  crack  are  identical. 

In  general  a  fatigue  crack  growth  law  can  be  formulated  by  coupling  the 
cyclic  stress-strain  distribution  ahead  of  the  crack  with  a  suitable  failure 
criterion. 

2.1  CRACK-TIP  FAILURE  CRITERIONS 

Researchers  in  the  past  have  used  different  failure  criterion  to  predict 
fatigue  crack  growth,  some  of  which  are  listed  below. 

2.1.1  Plastic  Strain  or  Total  Strain  Ahead  of  the  Crack 

In  this  approach  it  is  assumed  that  only  the  strains  above  the  yield  strain, 
Ey,  cause  the  damage  in  the  material.  Duggan  [9]  in  his  study,  has 
assumed  the  material  to  be  elastic  perfectly  plastic,  with  no  strain 
hardening  and  no  Bauschinger  effect.  The  failure  life  is  calculated  from 

a 

the  Coffin-Manson  [21]  law,  ACpN  f  =  C,  which  assumes  that  the  material 

will  fracture  when  the  cyclic  plasticity  builds  up  to  a  critical  value. 

Lehr  and  Liu  [10]  combine  the  Coffin-Manson  law  and  Miners  rule  [22] 
and  consider  a  volume  element  being  strained  with  increasing  amplitudes 
while  traversing  the  monotonic  plastic  zone,  rm. 
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With  a  strain  distribution 


where  Sy  is  the  yield  strain,  x  is  the  distance  from  the  crack  tip  and  p  is 

an  exponent.  Failure  is  assumed  to  occur  when  the  total  strain  exceeds  the 
fracture  strain. 

/ 

McClintock  [23]  assumes  that  for  a  crack  opening  in  Mode  III,  the 
average  plastic  shear  strain  in  the  volume  element  governs  the  crack 
propagation  process.  Failure  being  when  it  equals  the  plastic  fracture 
strain  for  Mode  I. 


2.1.2  Magnitude  of  the  Crack  Tin  Opening 

Some  researchers  have  used  the  crack  extension  per  load  cycle  equal  to 
the  crack-tip  opening  displacement,  Au,  caused  by  the  cyclic  stress  and 
hence  equal  to  the  crack  growth  rate. 

Lardner  [24]  applied  the  crack-tip  opening  model  developed  for  the  static 
load  by  Swinden  [27  in  24]  to  the  cyclic  loading.  But  this  model  does  not 
predict  the  true  shape  of  the  crack  growth  curve.  It  can  only  predict  the 
crack  growth  in  the  Paris  region,  with  some  fitting. 


2.1.3  Energy  Based  Failure  Criterions 
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In  the  energy  based  approaches  it  is  assumed  that  the  absorbed  energy 
during  cyclic  loading  causes  the  crack  initiation  and  propagation. 

Ellyin  and  Kujawski  [17,  18]  explain  that  during  cyclic  loading,  part  of 
the  mechanical  input  energy  is  rendered  irrecoverable.  A  major  part  of 
this  irrecoverable  energy  is  converted  into  heat,  and  the  remaining  part  is 
responsible  for  the  crack  initiation  and  growth.  This  energy  is  the  area  of 
the  hysteresis  loop  and  it  does  not  change  after  the  initial  few  cycles. 
They  calculated  this  absorbed  plastic  strain  energy  per  cycle,  as: 


AW  = 


1-n'^ 
1  +  n'y 


AaAs^ 


(3) 


where  n’  is  the  strain  hardening  exponent.  Act  is  the  stress  range,  Ae’’  is 
the  plastic  strain  range. 

Fatigue  damage,  in  general,  and  crack  propagation,  in  particular  can  be 
related  to  this  energy  term. 

They  express  the  stress  range  Ac  and  the  plastic  strain  range  Ag^,  in 


terms  of  number  of  reversals  to  failure  2Nf  by 

Ac  =  2aV  (2Nf  and  AgP  =  2eV  (2Nf (4) 

where  b  is  the  fatigue  strength  exponent,  c  is  the  fatigue  ductility 

•  1 
exponent,  Cf  is  the  fatigue  strength  coefficient  and  Sf  is  the  fatigue 

ductility  exponent.  Product  of  these  ranges,  a  measure  of  the  absorbed 
plastic  strain  energy  density  is  expressed  as: 
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AaAeP=4ffVe'f(2Nff''" 


(5) 


The  “process  zone”  where  the  stress  and  strain  ranges  have  a  very  low 
gradient  was  considered,  as  a  region  where  this  absorbed  plastic  strain 
energy  will  cause  failure  after  a  given  number  of  cycles. 

In  another  study  by  Pandey  and  Chand  [19],  the  criterion  of  specific 
energy  ‘Wc’  used  by  McEvily  [8  in  19],  has  been  modified  by  using  the 
Rice  superposition  method  for  cyclic  loading  to  modify  the  specific 
energy  term.  They  consider  the  entire  reversed  plastic  zone  as  the  area  for 
the  plastic  energy  dissipation,  with  a  premise  that  irreversible  plastic 
work  done  within  this  area  is  a  source  of  energy  dissipation.  A 
‘characteristic  length’  described  as  the  area  divided  by  the  circumferential 
length  of  the  cyclic  plastic  zone  is  used  to  get  the  unit  increment  for  crack 
propagation.  The  crack  advances  a  unit  increment  when  the  energy 
dissipation  in  this  region  equals  the  area  below  the  cyclic  stress-strain 
hysteresis  curve. 

A  cyclic-stress  strain  hysteresis  curve  is  shown  in  Fig. 8.  Point  ‘A’ 
corresponds  to  twice  of  the  cyclic  yield  stress  and  strain,  and  point  ‘B’  is 
the  point  corresponding  to  the  stress  and  strain  ranges  with  2Nf  =  1.  The 
material  behavior  is  described  by  Ramberg-Osgood  equation.  The  energy 
absorbed  till  failure,  Wc,  is  the  area  of  the  cyclic  stress-strain  range  curve. 
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We  =  Area  OABCO  =  Area  OFBEO  -  Area  OABFO  -  Area  CBEC 


Figure  8.  Cyclic  stress-strain  hysteresis  curve  showing  area  under 
curve  equivalent  to  the  energy  absorbed  till  failure,  Wc. 

It  has  been  noted  that  the  energy-based  failure  criterions  are  more  suitable 
than  others.  One  of  the  explanation  being  that  the  energy  based  criteria 
use  product  of  the  stress  and  plastic  strain  range,  while  in  other  criteria 
the  intrinsic  interplay  between  the  stress  and  strain  ranges  is  overlooked. 
Therefore,  the  criteria  that  utilize  both  the  stress  "and  plastic  strain  ranges 
give  better  results  than  only  stress  or  strain  criteria. 
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2.2  DAMAGE  ACCUMULATION 


From  the  phenomenological  point  of  view,  the  fatigue  crack  growth  has 
been  related  to  the  cyclic  fatigue  damage  process  of  material  elements,  at 
the  ‘process  zone’  ahead  of  the  crack  tip  [10-18].  Most  of  these 
approaches  consider  a  damage  element  near  the  crack  tip  which  is 
modeled  as  a  uniaxial  fatigue  specimen. 

Plastic  strains  are  the  cause  of  damage  accumulation  at  all  levels  of 
loading  [7,  8],  but  when  a  material  element  enters  the  reversed  plastic 
zone  it  is  not  subjected  to  fatigue  damage  immediately.  Fatigue  damage 
accumulates  only  when  the  plastic  strain  range,  Asp,  exceeds  a  threshold 
value,  Aso,  for  the  material.  As  soon  as  this  accumulated  damage  reaches  a 
critical  level,  the  material  element  fails. 

Cumulative  fatigue  damage  is  an  old,  but  not  yet  fully  resolved  problem. 
There  are  various  theories  to  calculate  the  accumulated  damage  under 
cyclic  loading.  These  theories  are  grouped  into  six  main  categories  by 
Fatemi  and  Yang  [25]  as: 

(a)  linear  damage  rules,  (b)  non-linear  damage  rule  and  two-stage 
linearization  approaches,  (c)  life  curve  modification  methods,  (d) 
approaches  based  on  crack  growth  concepts,  (e)  continuum  damage 
mechanics  models,  and  (f)  energy  based  theories. 

The  linear  damage  rule  cannot  account  for  the  load  sequence  and 
interaction  effects  due  to  its  linear  nature.  Though  many  damage  models 
have  been  developed,  unfortunately,  none  of  them  enjoys  universal 
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acceptance.  Each  damage  model  can  only  account  for  one  or  several 
phenomenological  factors,  such  as  load  dependence,  multiple  damage 
stages,  nonlinear  damage,  load  sequence  and  interaction  effects,  small 
amplitude  cycles  below  the  fatigue  limit,  and  the  mean  stress  effects. 

Due  to  the  complexity  of  the  problem,  none  of  the  predictive  models  can 
encompass  all  of  the  factors.  The  applicability  of  each  model  varies  from 
case  to  ease.  Consequently,  the  Palmgren-Miner  linear  damage  rule  [22]  is 
still  dominantly  used  in  design,  in  spite  of  its  major  shortcomings. 

Wu  and  Cotterell  [8]  in  their  approach  based  on  Dugdale  model  and 
plastic  strain  energy  failure  criterion  similar  to  Rice,  inodel  the  damage 
accumulation  process  according  to  the  Miner’s  linear  damage  rule. 

Duggan  [9]  defines  the  damage  zone  ahead  of  the  crack,  Rd,  as  the  radius 
from  the  crack  tip  to  the  point  corresponding  to  the  endurance  strain  range 
He  also  assumes  that  the  damage  accumulates  linearly  according  to 
Miner’s  law. 

Glinka  [10-12]  modeled  the  fatigue  crack  growth  for  steel  and  aluminum 
alloys  by  using  local  stresses  and  strains  at  the  crack  tip,  and  the  damage 
accumulation  process  in  the  highly  strained  region  ‘process  zone’  in  the 
vicinity  of  the  crack  tip. 

In  one  of  the  most  recent  of  approaches  based  on  the  critical  plastic  strain 
energy  dissipation  by  Pandey  and  Chand  [19],  damage  accumulation  ahead 
of  the  crack  tip  has  been  successfully  modeled  by  using  the  Miner’s  law. 
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The  results  reported  in  the  above  studies  suggest  that  the  Palmgren- 
Miner’s  linear  damage  rule  is  the  most  popular  damage  accumulation 
approach,  which  can  be  easily  used  in  conjunction  with  variety  of  failure 
criteria. 
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2.3  EFFECT  OF  APPLIED  R-RATIO  ON  FAIGUE  CRACK  GROWTH 
The  load  ratio  or  R-ratio,  is  defined  as  the  ratio  of  Smin/S^ax  or  Kmin/K^ax, 
where  the  S  is  the  stress  and  K  is  the  stress  intensity  factor. 

It  is  known  that  an  increase  in  the  R-ratio  for  the  cyclic  loading  at  a  given 
AK  causes  the  crack  growth  rate  to  be  higher.  This  R-ratio  effect  on 
fatigue  crack  growth  is  shown  in  Fig. 9. 

For  some  metals  the  R-ratio  effect  is  small  or  is  non  existent  in  the  Paris 
region.  The  R-ratio  effects  are  more  pronounced  in  the  near  threshold 
region  and  in  the  rapid  unstable  crack  growth  region. 
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Figure  9.  Long  crack  fatigue  data  for  2024-T351  Al.  alloy  [26]. 
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A  number  of  investigators  have  studied  the  R-ratio  effects  on  near 
threshold  region  in  laboratory  air  and  vacuum,  in  terms  of  relation 
between  AKth  and  R.  The  general  trend  shows  that  for  many  materials  in 
air,  AKth  tends  to  decrease  with  increasing  R,  but  for  some  materials  and 
R  >  0.5,  the  sensitivity  is  less.  Most  experimental  results  in  vacuum 
indicate  that  R-ratio  effect  on  the  AKth  is  less  pronounced  or  diminished 
[3]. 
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2.4  CONCEPTS  FOR  CORRELATION  OF  R-RATIO  EFFECTS 


There  are  several  concepts  that  have  been  adopted  to  account  for  the  R- 
ratio  effect,  namely  a)  Crack  closure,  b)  Residual  compressive  stresses  c) 
Environmental  influence.  All  the  approaches  tend  to  calculate  an 

‘effective  stress’,  a ,  or  ‘effective  stress  intensity  factor’,  K ,  to  account 
for  the  R-ratio  effect  on  the  fatigue  crack  growth. 


In  1970  Walker  [26]  used  an  effective  stress  approach 


a 


m 


(6) 


which  was  able  to  correlate  effects  of  R-ratio  on  the  fatigue  life  for  7075- 
T6  and  2024-T3  A1  alloys. 


For  the  case  of  fatigue  crack  growth,  the  Eq.(6)  was  modified  to  the 
following  form 


■  K  =  Kl-)AK"’  (7) 

which  correlated  the  crack  growth  data  for  positive  R-ratios  fairly  well.  In 
his  approach  ‘m’  was  assumed  to  be  a  material  parameter  which  was 
calculated  from  experimental  data. 

In  1971  Elber  [27]  introduced  a  crack  closure  concept,  in  which  he  stated 
that  the  crack  is  open  during  only  some  part  of  the  load  cycle.  He  defined 
Kop  as  the  stress  intensity  factor  at  which  the  crack  opens,  while  cycling 
from  Kmin  to  Kmax-  It  is  assumed  that  there  is  no  damage  ahead  of  the 
crack  tip  for  the  portion  of  the  cycle  in  which  crack  is  closed.  Therefore, 
the  effective  stress  intensity  factor  range,  Keff  is  calculated  as: 
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^^eff  ^max  ^op  (8) 

Fig.  10  illustrates  the  typical  relation  among  K^ax,  K^in  and  Kop  for  load 
cycles  with  R  =  0.1  and  0.5. 


Figure  10.  Load  cycle  diagram  for  R=0.1  and  R=0.5. 

The  corresponding  Keff  is  calculated  as: 

For  R=0.1  ;  Kop  >  Kn,in  ;  AKeff  =  K„ax  -  Kop  =  7 
For  R  — 0.5  ,  Kop  <  Kmln  i  AKeff  ~  K^ax  ~  Kmin  “5 

Fiber  suggested  that  residual  tensile  displacements  left  in  the  wake  of  the 
crack  tip  plastic  zone  cause  the  crack  to  remain  closed  even  when  the 
applied  load  is  positive. 

Since  then  the  concept  of  crack  closure  has  been  most  widely  adopted  as  a 
critical  mechanism  responsible  for  the  R-ratio  effects. 

Dinda  and  Kujawski  [28]  used  a  fatigue  crack  driving  force  parameter 
similar  to  Walker’s  approach,  difference  being  that  it  is  calculated  solely 
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using  the  positive  part  of  the  range  of  applied  stress  intensity  factor,  AK'^, 
and  the  corresponding  maximum  value  of  SIF,  Kmax- 

aKd  =  (9) 

Using  this  approach  the  R-ratio  effects  on  fatigue  crack  growth  were 
successfully  correlated  for  a  number  of  materials  [28]. 

Another  approach  by  Harter  [29],  called  the  Harter-T  method  is  an 
adaptation  of  the  Walker  equation  taken  on  point-by-point  basis.  In  this 
implementation  25  da/dN  and  corresponding  AK  values  are  used  to  define 
the  crack  growth  behavior.  The  tabular  data  utilizes  the  Walker  equation 
on  a  point-by-point  basis  (Harter  T-Method)  to  extrapolate/interpolate 
data  for  any  R  value.  Illustration  in  Fig.  11  represents  the  tabular  data  for 
two  R-ratios  which  are  used  in  Walker  equation,  Eq.(7),  to  determine  m. 


Figurel  1.  Illustration  of  point  by  point  basis  Harter-T  method. 
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Since  da/dN  -  C[AK(1-R)^™''^]"  (  Walker  equation) 

At  a  given  da/dN,  data  at  any  two  positive  R-values  becomes: 

AK,(1-Ri)<”’-')=AK2(1-R2)<'"“"  (10a) 

f  AKi  ^  /  f  1“R 

m  =  l+  loglO  /loglO— ^  (10b) 

When  prediction  is  made  for  the  negative  R-ratios,  a  distinction  is  made 
such  that  Kmax  is  used  in  place  of  AK. 

According  to  Vasudevan  and  Sadananda  [30]  the  fatigue  crack  growth  is 
based  on  two  basic  fracture  mechanics  parameters  Kmax,  AK  and  the 
internal  stress  contribution  to  Kmax-  They  propose  that  the  internal 
stresses  are  the  missing  link  that  can  bridge  the  gap  between  four  main 
stages  of  damage  from  nucleation  to  final  failure. 
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2.5  DRAWBACKS  OF  THE  CRACK  CLOSURE  METHODOLOGY 


After  Elber  [27]  proposed  the  crack  closure  model,  much  focus  was 
towards  relating  the  load  interaction  effects  with  the  plasticity  induced 
crack  closure  operating  in  the  crack  wake.  Over  the  time  many  researches 
have  shown  that  other  types  of  crack  closure  such  as  a)  surface  roughness 
or  asperity  induced  closure,  b)  fracture  surface  mismatch  induced  closure 
c)  fracture  surface  oxidation  induced  closure  and  d)  debris  related  closure 
have  to  be  considered. 

In  1993  Sadahanda  and  Vasudevan  [31]  reconsidered  the  plasticity- 
induced  closure  and  its  effects  on  FCG  behavior.  They  argued  that,  since 
the  closure  occurs  behind  the  crack  tip,  it  has  rather  limited  effect  on 
damage  process  that  takes  place  in  front  of  the  crack.  They  demonstrated 
that  there  is  no  significant  contribution  to  closure  due  to  plastic  stretch  in 
the  crack  wake.  Further,  when  asperity  or  roughness  induced  closure  is 
present,  the  contribution  is  small,  approximately  one  quarter  of  that 
determined  from  compliance  measurements. 

Recently  in  2003  Silva  [32]  showed  that  the  closure  concept  is  inadequate 
at  the  negative  stress  ratios.  He  conducted  experiments  at  negative  stress 
ratios  to  focus  on  the  crack  tip  plasticity  and  surface  roughness  effects  on 
closure.  Study  confirmed  that,  at  negative  stress  ratios,  the  closure 
changes  with  Pmax  for  the  same  R-ratio.  He  concluded  that  the  damage 
accumulation  principles  or  residual  stress  concepts  seem  to  be  more 
adequate  to  explain  R-ratio  effects  than  the  crack  closure  concept. 
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Original  crack  closure  concept  assumes  that  the  whole  crack  opens 
instantly  at  opening  load,  Pop.  Thus,  at  any  load  P<Pop  crack  is  closed  and 
fully  shielded  from  the  damage. 

But,  for  most  metals,  crack  growth  in  the  near  threshold  region  is 
associated  with  a  single  shear  mode  of  crack  growth,  which  gives  rise  to 
faceted  fracture  surfaces  and  their  associated  mismatch.  Therefore,  in  near 
threshold  regions  the  roughness  or  mismatch  dominated  closure  would  be 
expected. 

Donald  and  Paris  [1)3,  34]  modeled  this,  using  so  called  partial  crack 
closure  model,  where  closure  of  crack  faces  only  partially  shields  the 
crack  tip  from  fatigue  damage.  They  proposed  the  following  modification 
to  the 

AKeff  relation  in  case  of  the  partial  closure  in  near  threshold  region 

AK™=K„„--K^  if  -K„p>K„,„  (Ila) 

^^2/?!  =K^ax“^^min  ~^op-K^min  (Hb) 

Usually  Pop  is  determined  according  to  the  ASTM  standard  E647 
recommendations,  which  corresponds  to  the  load  that  causes  the  2% 
deviation  in  the  slope  of  a  load-displacement  curve.  Then  the  Keff  is 
calculated  from  this  Pop  value.  However,  it  is  often  observed  that  such 
calculated  Keff  for  a  low  stress  ratio  (e.g.  R=0.1)  at  threshold  is  lower 
than  the  applied  AK  for  a  high  stress  ratio  (R>0.7  or  higher)  where  crack 
closure  is  usually  absent. 
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This  contradiction  is  illustrated  by  an  example  below  [41] 

Fig.  12  shows  a  typical  data  for  crack  growth  rate  (da/dN  vs.  AK)  at  two 
stress  ratios  R=0.7  and  0.1. 


Figure  12.  Illustration  of  typical  crack  growth  data  for  R=0.7 
and  R=0. 1 . 

Assuming  a  similar  kind  of  loading  as  shown  in  Fig.  10  with  Kop=6 
At  load  ratio  R=0.1;  K^ax  =10  and  K^in  =1, 

At  load  ratio  R=0.7;  Kmax  =10  and  Kmin  =7, 

and  the  opening  stress  intensity  factor  Kop  =  6 

Therefore  using  Eqn.(8): 

for  R^O.l  ,  Kop  ^  Kmin  >  AKeff  =  Kinax"Kop  “4 

for  R— 0.7  5  Kop  <  Kmin  >  AKeff  =  Kmax'Kmin  =  3 

So  theoretically  the  data  should  correlate  as  shown  in  the  Fig.  13. 
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Figure  13.  Expected  theoretical  correlation  using  the  Keff  method. 

However  it  has  been  shown  by  many  authors  that  using  the  “crack  closure 
methodology”  the  correlation  shown  in  Fig.  14  is  obtained,  where  the 
curves  for  R=0.1  and  0.7  are  in  reversed  order  than  in  Fig.l3. 

This  reverse  order  of  curves  indicated  that  the  effective  threshold  for 
R=0.1  is  smaller  than  for  R=0.7.  This  contradicts  the  crack  closure 
concept  because,  the  data  for  R=0.1  with  crack  closure  has  lower 
threshold  than  the  data  for  R=0.7  where  closure  is  absent. 


Figure  14.  Correlation  obtained  using  crack  closure  methodology. 
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The  above  mentioned  contradiction  has  been  observed  by  many  authors 
and  Donald  and  Paris  [33,  34]  have  explained  it  on  the  basis  of  the  partial 
crack  closure  phenomenon  discussed  earlier. 

One  of  the  most  important  drawback  of  the  crack  closure  methodology  is 
the  difficulty  of  measuring  the  crack  opening  load.  The  crack  opening 
load  depends  on  the  measurement  location  and  the  technique  used.  Fig.  15 
shows  a  typical  load  vs.  displacement  curve  which  indicates  that  there  is  a 
gradual  transition  from  fully  closed  to  fully  open  crack. 

But  the  closure  method  takes  a  single  point  as  Pop  below  which  the  crack 
is  fully  closed  and  above  it  is  fully  open.  To  determine  this  point  there  are 
several  different  methods  that  can  be  applied,  and  each  gives  a  slightly 
different  result  [41]. 


Figure  15.  Load  vs.  displacement  curve  showing  compliance  change. 
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Me  Clung  [35]  after  doing  extensive  review  on  crack  closure  concluded 
that  there  are  three  distinct  regimes  of  crack  closure  namely  a)  near 
threshold,  b)  Paris  region  and  c)  high  crack  growth  region. 

Donald  [36]  also  showed  that  the  fatigue  crack  growth  rate  is  not 
determined  solely  by  AKeff  but  also  depends  on  the  Kmax- 
With  all  the  drawbacks  about  the  crack  closure  and  the  conclusions 
mentioned  by  different  authors  it  can  be  inferred  that  there  is  a  need  to 
develop  a  fatigue  crack  growth  model  that  has  both  the  Kmax  and  AK 
effects  on  the  crack  growth  rate.  The  model  should  take  into  account  not 
only  what  is  happening  behind  the  crack  tip,  but  also  what  is  happening 
ahead  of  the  crack  tip. 
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2.6  TWO  PARAMETER  (Kmax  and  AK)  UNIFIED  APPROACH 
In  the  work  done  by  Vasudevan  and  Sadananda  [3.7-39],  they  explain  that 
there  are  two  parameters,  AK  and  K^ax,  that  are  required  to  describe  the 
crack  driving  force  that  controls  the  fatigue  crack  growth.  This  approach 
has  been  given  the  name  of  ‘unified  approach’  to  fatigue  crack  growth. 

In  early  approaches  starting  from  Paris  et  al.  [6],  AK  and  R  have  been  the 
two  parameters  to  quantify  the  fatigue  crack  growth  rate.  The  existence  of 
threshold  in  terms  of  AK  has  been  considered,  as  it  is  a  contributing  factor 
to  the  driving  force.  Threshold  for  R  is  not  considered,  since  it  is  not  a 
driving  force.  Hence,  the  effects  of  R-ratio  on  fatigue  crack  growth  have 
been  explained  mainly  by  the  crack  closure  approach. 

In  contrast,  the  unified  approach  considers  Kmax  also  to  be  a  required 
parameter  to  quantify  the  driving  force  in  addition  to  AK.  The  existence 
of  the  Kmax  threshold  has  been  reported  by  the  authors,  which  is  a  definite 
indication  that  Kmax  is. a  crack  tip  driving  force. 

Vasudevan  et  al.  [30]  have  also  shown  that  all  deviations  from  steady 
state  long  crack  growth  behavior  can  be  accounted  for  by  presence  of 
“internal  stress”  gradients  that  introduce  internal  stress  intensity  factor 
Kint-  The  Kint  mainly  affects  the  Kmax  and  causes  perturbations  in  it.  The 
value  of  AK  is  not  much  affected  since  Kjnt  acts  on  both  Kmax  and  Kmin- 
The  unified  approach  suggests  the  following  procedure  for  fatigue  crack 
growth  modeling: 

(a)  Analytical  models  must  include  Kmax,  AK  and  Kint. 
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(b)  Understanding  and  estimating  the  internal  stresses  is  fundamental  to 
the  development  of  a  reliable  life  prediction  model. 

(c)  Modeling  must  account  for  the  accumulation  of  damage  from  crack 
nucleation  stage  to  short  crack  to  long  crack  to  final  failure. 

(d)  Physics-based  descriptions  of  fatigue  crack  nucleation  and  growth  are 
needed  in  order  to  obtain  a  reliable  life  prediction  model. 
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3.  SCOPE  OF  THE  WORK 

The  purpose  of  this  study  is  to  develop  a  mathematical  model  for  the 
fatigue  crack  propagation  utilizing  the  stress  and  strain  state  near  the 
crack  tip  region  and  damage  accumulation  process.  The  methodology  will 
be  based  on  the  use  of  both  Kmax  and  AK  as  the  contributing  parameters  to 
the  driving  force  and  to  understand,  how  these  parameters  affect  the 
damage  accumulation. 

Effort  will  also  be  made  to  understand  and  estimate  the  internal  stresses, 
and  to  model  the  effect  of  these  internal  stresses  and  Kjnt  on  the  applied 
Kmax  • 

This  will  not  only  help  to  predict  the  crack  propagation  rates  without  the 
experimental  determination  of  empirical  parameters,  but  it  would  also 
lead  to  a  better  understanding  of  the  effect  of  these  factors  on  the  fatigue 
crack  growth  behavior. 

The  crack  tip  stress-strain  distributions  within  and  beyond  the  plastic 
zone  will  be  evaluated  adopting  Rice  [7]  and  the  elastic  solution, 
respectively.  It  has  been  widely  accepted  [8-12,  15-18,  20]  that  there  is  a 
small  region  ahead  of  the  crack,  where  due  to  yielding  and  crack  tip 
blunting  the  stress  and  strain  magnitudes  are  finite,  so  the  singularity  at 
crack  tip  inherent  to  these  models  can  be  avoided.  This  small  region  in 
front  of  the  crack  tip  is  termed  the  ‘process  zone’,  p*.  Subsequently,  the 
damage  accumulation  in  the  process  zone  will  be  determined  and  analyzed, 
in  relation  to  how  ifgoverns  the  progressive  growth  of  the  crack.  In  this 
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study  it  will  be  assumed  that  the  damage  accumulates  according  to  the 
Palmgren-Miner  [15]  linear  damage  rule. 

The  effect  of  the  internal  stresses  on  K^ax  will  be  modeled  by  evaluating 
Kint  using  the  weight  function  methodology  proposed  by  Glinka  [54-57], 

A  graphic  user  interface  (GUI)  will  also  be  developed  based  on  the  model, 
which  will  be  programmed  in  MATLAB.  Finally,  the  model  will  be 
verified  by  comparing  the  predictions  of  the  model  to  the  data  available  in 
literature  for  a  spectrum  of  materials.  The  verification  of  the  predictions 
by  the  current  model  will  consist  of  following  steps: 

(a)  Predict  the  crack  growth  data  for  the  constant  amplitude  loading  and 
compare  it  to  the  available  data  published  in  literature. 

(b)  Next,  the  predictions  of  the  weight  function  subroutine,  developed  in 
MATLAB,  will  be  compared  to  the  predictions  by  Glinka’s  FALPR 
application. 

(c)  Internal  stresses  ahead  of  the  crack  tip  for  the  overload  will  be 
utilized  to  predict  the  Kjnt  profile. 

(d)  The  predictions  for  the  overload  effected  crack  growth  data,  by  Kim 
modeling,  will  be  compared  to  the  data  available  in  the  literature. 

(e)  Predictions  for  the  R-ratio  effect  by  the  current  model  will  be 
compared  to  the  data  available  in  the  literature. 
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4.  STRESS  AND  STRAIN  DISTRIBUTION  AHEAD  OF  THE  CRACK  TIP 


Rice  [7]  analytically  solved  the  stress  and  strain  distribution  ahead  of  the 
crack  tip  for  a  stationary  crack  in  anti-plane  shear  (mode  III  loading) 
under  small  scale  yielding.  The  corresponding  shear  stress  and  shear 
strain  relationship  was  assumed: 


t  =  G7  for  7<7o,T<To 


And 


T=Tr 


f 
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for  7>7o,T>To 


(18) 


Where  Tq  and  7^  are  the  yield  strength  and  strain  respectively  and  n  is  the 
strain  hardening  exponent. 

The  following  equations  for.  mode  III  were  derived  [7] 
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A  similar  analysis  for  tensile  loading  (mode  I)  is  not  available.  However, 
McClintock  [23]  had  suggested  that  there  is  the  analogy  between  mode  III 
and  mode  I  for  a  case  where  displacements  parallel  to  the  crack  are  small 
as  compared  to  those  normal  to  the  crack  surface. 
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Assuming  the  following  stress  and  strain  relationship 
a  =  Es  for  s  <  Sq,<j  < 

And 


(20) 
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The  analogous  equations  for  Mode  I  are 


(T  =  CT„ 

max  o 


■^max 

2 

(1  +  n);rcrQX 

> 

£  *“  £ « 

max  0 

2 

(1  +  n);rcrQX 

q+n 


(21) 


Assuming  that  the  plastic  strain  components  at  each  point  within  the 
plastic  zone  remain  proportional  to  each  other  the  following  equations  can 
be  adopted  for  unloading 

^1/ 


Act  =  2cr^ 
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The  product  of  stress  and  strain  ranges  is  given  by  multiplying  the  above 
two  equations 

■  2 
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The  cyclic  plastic  zone  size,  rc,  can  be  approximated  by  setting 
Act  =  2<Jq,As  ^ISq  in  Eq.(23). 
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Similarly,  by  setting  <7  =  (To  Eq.(23)  the  monotonic  plastic  zone,  ‘Tp’ 
can  be  approximated  as 


(l  +  n)7tal 


It  has  been  shown  by  Kujawski  and  Ellyin  [18]  that  Eq.(21)  and  Eq.(22) 
are  only  valid  inside  the  monotonic  and  cyclic  plastic  zones  respectively, 
whereas  the  elastic  solution  outside  the  plastic  zones  should  be  used. 
Assuming  plane  stress,  the  elastic  stress  and  strain  and  their  ranges  can  be 
calculated  from  the  elastic  solution  given  by: 


(26) 


AK  AK 
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While  calculating  the  elastic  solution  from  the  Eq.(26)  and  Eq.(27),  the 
distance  ‘x’  should  be  offset  by  approximately  half  of  the  monotonic 
plastic  zone  size  for  the  loading  part,  Om,  and  by  approximately  half  of  the 
cyclic  plastic  zone  for  the  unloading  part,  Oc.  This  offset  is  applied  to 
make  the  prediction  to  be  continuous  at  the  monotonic  and  cyclic  plastic 
zones. 


The  exact  values  of  these  offsets  can  be  derived  as  follows: 

At  X  =  tp,  the  value  of  stress  by  both  the  rice  solution  and  elastic  solution 
should  be  equal  to  the  yield  stress,  Gq. 

Substituting  CTmax  =  CTo  and  solving  Eq.(21)  for  x. 
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Once  both  the  loading  and  the  unloading  profiles  have  been  calculated, 
the  stress  profile  for  the  minimum  load  is  found  out  by  subtracting  the 
solution  for  unloading  from  the  solution  for  loading.  In  this  study,  the 
Rice  solution  within  the  plastic  zone  and  the  elastic  solution  beyond  this 
zone  were  adopted  for  the  stress  and  strain  calculation  in  the  direction 


perpendicular  to  the  crack  plane. 

Since  both  solutions  exhibit  a  singularity  as  x->  0,  there  exists  a  small 
region  ahead  of  the  crack  tip  called  ‘process  zone’  where,  due  to  yielding 
and  crack  tip  blunting,  the  stresses  and  strains  have  a  finite  magnitude. 
The  stress  and  strain  gradient  in  reality  is  much  smaller  in  this  region 


than  that  predicted  by  analytical  solutions.  So  for  the  predictions  to  be 
very  close  to  the  real  behavior,  the  stresses  and  strains  are  assumed  to  be 
constant  in  the  process  zone.  Fig.  16  illustrates  the  stress  and  strain 
distribution  ahead  of  the  crack  in  mode  I  loading  for  4140  steel  for  R  =0.1 
at  two  different  load  levels  and  Fig.  17  illustrates  the  stress  and  strain 
distribution  ahead  of  the  crack  tip  for  4140  steel  for  R=0.5.  Both  the 
graphs  have  been  plotted  for  a  distance  ahead  of  the  crack  tip,  equal  to 
twice  the  monotonic  plastic  zone  size.  As  it  was  described  above,  the 
stress-strain  distributions  within  and  beyond  the  plastic  zone  have  been 
evaluated  adopting  Rice  and  the  elastic  solution,  respectively. 


Figure  16.  Stress  distribution  ahead  of  the  crack  tip  for  the  4140 
steel  at  two  applied  AK  values  for  R=0.1. 
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Stress  fields  for  Rice  +  Elastic  solution 


Figure  17.  Stress  distribution  ahead  of  the  crack  tip  for  the  4140 
steel  at  two  applied  AK  values  for  R=0.5. 


As  it  can  be  seen  from  the  Fig.  16  and  Fig. 17,  the  change  in  R-ratio  from 
0.1  to  0.5  with  same  AK,  changes  only  the  loading  part  of  the  profile  and 
the  monotonic  plastic  zone  size.  The  unloading  part  of  the  profile  and  the 
cyclic  plastic  zone  sizes  are  the  same  in  both  cases. 

Also  to  be  noted  is  the  fact  that  very  far  away  from  the  crack  tip  the  local 
R-ratio  reaches  the  applied  R  ratio.  The  corresponding  mechanical 
properties  for  4140  steel  are  listed  in  Table  1. 
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Table  1.  Mechanical  properties  of  4140  and  4340  steel  [15,  46,  47] 


PROPERTY 

SYMBOL 

0050 A  steel 

4340  STEEL 

UNITS 

Young's  Modulus 

E 

205 

209 

GPa 

Yield  Stress 

(TO 

645 

724 

MPa 

Cyclic  hardening  exponent 
Cyclic  strength  coefficient 

n* 

H' 

0.15 

1640 

0.146 

1720 

MPa 

Fatigue  strength  coefficient 

Of 

1530 

1713 

MPa 

Fatigue  ductility  coefficient 

Cf 

0.65 

0.83 

Fatigue  strength  exponent 

b 

-0.087 

-0.095 

Fatigue  ductility  exponent 

c 

-0.6 

-0.65 

Table  2.  Mechanical  properties  of  2324-T39  and  7475  Al.  alloy  [59]. 


PROPERTY 

SYMBOL 

2324-T39  Al. 

7475  Al. 

UNITS 

Young's  Modulus 

E 

71. 

71 

GPa 

Yield  Stress 

(TO 

430 

460 

MPa 

Cyclic  hardening  exponent 

n' 

0.106 

0.11 

Cyclic  strength  coefficient 

H' 

877 

977 

MPa 

Fatigue  strength  coefficient 

_ o-f 

927 

1466 

MPa 

Fatigue  ductility  coefficient 

ef 

0.4 

0.262 

Fatigue  strength  e7q)onent 

b 

-0.0642 

-0.4985 

Fatigue  ductility  exponent 

c 

-0.4009 

-0.088 
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5.  FAILURE  CRITERION  FOR  SMOOTH  SPECIMEN 


Most  -engineering  components  and  structures  experience  both  alternating 
stresses  ‘Oa’  and  a  mean  stress  ‘Om’  during  service  life.  Often  defects  or 
crack-like  flaws,  ranging  from  micro  to  macro  sizes  may  be  present  which 
cause  the  stress-strain  concentration.  The  local  stress-strain  field  in  the 
neighborhood  of  the  crack  governs  the  damage  accumulation  [7-13].  To 
correlate  fatigue  life,  a  proper  damage  parameter  is  needed  that  has  an 
appropriate  measure  of  mean  stress  effect  and  the  applied  stress-strain 
amplitude.  The  SWT  parameter  or  the  Morrow  Equation  can  be  used  as  a 
failure  criterion  to  predict  the  failure  life  of  the  smooth  specimen.  In  this 
analysis,  the  SWT  parameter  has  been  adopted  to  calculate  the  failure  life 
(Nf)  of  each  element  ahead  of  the  crack  tip.  The  whole  specimen  is 
divided  into  the  small  elements  of  width  p*  (see  Fig. 18)  and  each  element 
.  is  considered  as  a  smooth  uniaxial  specimen  for  which  the  fatigue  life  can 
be  estimated  using  the  SWT  parameter  or  Morrow  equation.  The  crack 
advances  by  one  element  when  the  failure  criterion  for  that  element  is  met. 
The  failure  criterion  is  based  on  damage  accumulation,  and  is  explained  in 
the  sections  ahead. 
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Figure  18.  Representation  showing  the  elements  ahead  of  the  crack 
tip  considered  as  smooth  specimens. 


5.1  SMITH-WATSON-TOPPER  (SWT)  PARAMETER 
The  SWT  equation  for  a  smooth  specimen  can  be  written  as 

.,2 


Af  CT  f  OK 

^max^a  ~  ^max  «  “  (2Nf )  +  cr'f  £'j-'(2Njp) 

^  E 


b+c 


(32) 


Since  a. 


Act 


max 


2(1  -  R) 


,  by  substituting  this  into  Eq.(32) 


In  a  modified  form  SWT  can  be  written  in  terms  of  the  stress  and  strain 
ranges  as 


AcjA£‘ 


..2 


2(1-R^*)  E 


=  -^(2Nf)2‘’  +  c7Vff'(2Nf)^-'‘^ 


(33) 
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Eq.(33)  was  obtained  using  =7—^  where  Rp  =  (o-^n  / ^-^ax L > 

A 

the  stress  ratio  in  the  process  zone,  p*. 

Eq.(33)  results  in  a  single  curve,  which  can  be  drawn  for  various  values  of 
Nf  vs.  the  SWT  parameter,  LHS  ofEq.(33). 

The  failure  life  Nf  can  be  evaluated  faster  by  interpolating  the  values  from 
the  SWT  curve  rather  than  solving  the  equation  numerically. 

The  RHS  of  Eq.(33)  has  two  parts: 

The  elastic  part,  and  the  plastic  part,  cTf  s\ 

E 

Fig. 19  shows  the  modified  SWT  curve  for  4140  steel  consisting  of  250 
values  of  Nf  equally  spaced  on  a  log-log  plot. 


Figure  19.  The  plot  of  the  modified  SWT  curve  for  4140  steel. 
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Slope  of  the  elastic  part  and  plastic  part  is  different  for  different 
materials,  and  depend  on  the  exponents  ‘b’  and  ‘c’.  Also  the  transition 
point  between  the  elastic  part  and  plastic  part  is  important  in  the  crack 
growth  prediction,  which  will  be  discussed  in  detail  in  the  sections  ahead. 


Table  3.  Mechanical  properties  of  A533-B1  and  C-Mn  steel  [15]. 


SYMBOL  A533-B1  STEEL  C-Mn  STEEL  j  UN 


Young's  Modulus  _ _ _ E _ [ 

Yield  Stress 


PROPERTY 


Cyclic  hardening  e^onent 
Fatigue  strength  coefficient 
Fatigue  ductility  coefficient 
Fatigue  strength  emonent 


Fatigue  ductility  e:iq)onent 


C-Mn  STEEL  i 

UNITS 

208  1 

GPa 

372  1 

MPa 

0.141 

868  i 

1  MPa 

0.15 

-0.514 

-0.101* 

5.2  MORROW  EQUATION 

The  approach  suggested  by  Morrow  can  be  expressed  as  an  equation 
giving  the  equivalent  completely  reversed  stress  amplitude,  Oar,  which  is 
expected  to  produce  the  same  life  as  given  by  combination  of  stress 
amplitude  Oa  and  mean  stress  Om. 


<^01  = 


(34) 


1 

^  » 


The  stress  life  equation  for  completely  reversed  loading  is 

C7„  =  f;(2Nf)^ 


(35) 


By  combining  Eq.(34)  with  Eq.(35),  where  cjf  and  b  are  evaluated  for  zero 
mean  stress  or  R-ratio  =  - 1 . 


C7a  =0-f 


1 

r  \z 

1 

V 


"|b 


(2N,) 


(36) 


Comparing  it  to  Eq.(35)  the  effect  of  mean  stress  on  Nf  is  such  that 

r 


1 


N  =Nf 


1 


V 

Applying  same  modification  to  the  strain-life  equation 


(37) 


=a(2N-)’’+.;(2N*r 


(38) 
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By  substituting  N*  value  into  Eq.(38)  and  neglecting  the  mean  stress 
effect  on  the  plastic  strain  part  (second  term)  in  Eq.(38)  a  single  equation 
for  family  of  strain-life  curves  can  be  obtained. 

+  4(2Nf  )'  (39) 

The  convenient  graphical  procedure  for  solving  this  equation  can  be, 
plotting  the  strain-life  curve  at  zero  mean  stress.  And  for  any  mean  stress, 
first  the  N*  values  can  be  read  from  the  curve,  which  can  be  solved  with 
Eq.(37)  to  get  the  Nf  value  for  any  mean  stress. 

This  kind  of  approach  with  graphical  procedure  is  shown  in  the  Fig. 21 
below: 


01 


o 

I 

Q. 

<  0,01 
C 
1 
to 


0.001 

1.E+00  l.E-K)1  1.E-K)2  1.E+03  1.E+04  1,E-K)5  1.E+06  1.E+07  1.E+08 

N*  Cycles  to  failure 

Figure  21.  Mean  stress  data  for  AISI  4340  Steel,  plotted  vs.  N* 
according  to  the  Morrow  equation. 
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If  the  family  of  strain-life  curves  is  plotted  it  would  look  like  shown  in 
Fig. 22,  but  it  will  be  a  long  and  time  consuming  method  to  interpolate  the 
values  between  this  family  of  the  curves. 


Figure  22.  Family  of  strain  life  curves  for  4340  Steel  plotted  vs.  Nf 
according  to  the  Morrow  equation. 


Both  the  approaches  SWT  and  Morrow  equation  are  in  current  use,  and  no 
consensus  exists  that  any  one  of  them  is  superior  to  the  other.  The  data 
shown  in  Fig. 21  is  plotted  in  Fig. 23  using  the  SWT  equation. 
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1.E+00  1.E+01  1.E+02  1.E+03  1.E+04  1,E+05  1.E+06  1.E+07 

Nf,  Cycles  to  failure 

Figure  23.  Mean  stress  data  for  AISI  4340  Steel,  plotted  vs.  Nf 
according  to  the  SWT  parameter. 

Morrow  approach  gives  good  results  for  steels,  and  in  some  cases  it  gives 
better  results  than  the  SWT  parameter.  However  for  some  aluminum  alloys 
SWT  gives  better  results. 

Dowling  [3]  concludes  that  SWT  parameter  appears  to  give  good  results 
for  wide  range  of  materials  and  is  a  good  choice  for  general  use.  However, 
it  may  give  non-conservative  estimates  for  compressive  mean  stresses. 
Where  refinement  is  desired,  the  choice  among  these  parameters  can  only 
be  based  on  fatigue  test  data  involving  mean  stresses.  In  the  present  study 
SWT  parameter  will  be  used  for  all  the  calculations  and  driving  force 
derivations. 
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6.  RELATION  BETWEEN  SWT  AND  AK*  PARAMETER 


Recently,  Glinka  [40]  derived  the  relation  between  the  exponent  a  in  the 

two  parameter  crack  driving  force  equation  AK  =  (Kjjj3x)“(AK 

proposed  by  Kujawski  [41]  and  the  cyclic  strain  hardening  exponent  n’ 
utilizing  the  SWT  parameter.  He  used  the  generalized  Neuber  law  for  the 
stress  and  strain  calculation  in  the  material  elements  at  the  crack  tip 
region.  A  similar  analysis  is  presented  in  this  section  where  the  stresses 
and  strains  in  the  crack  tip  region  are  calculated  using  modified  Rice  and 
elastic  solutions  discussed  in  the  previous  section. 

The  un-cracked  region  of  the  specimen  is  divided  into  small  elements  of 
width  equal  to  p*  as  shown  in  Fig. 24.  It  is  assumed  that  due  to  cyclic 
loading  the  significant  fatigue  damage  is  accumulated  predominately  in 
the  first  element  of  width  p*  called  the  process  zone.  When  the 
accumulated  fatigue  damage  in  the  process  zone  reaches  a  critical  value 
according  to  the  SWT  parameter,  the  crack  tip  will  advance  to  the  next 
element  and  the  process  is  repeated.  In  other  words,  it  is  assumed  that  the 
material  in  the  process  zone  can  be  considered  as  a  hypothetical  uniaxial 
specimen. 
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Figure  24.  Schematic  of  the  specimen  divided  into  elements 
of  width  p*. 


It  can  be  assumed  that  in  the  Paris  region  of  crack  growth  the  local 
deformation  in  the  crack  tip  process  zone  is  dominated  by  the  plastic 
strain. 

The  corresponding  plastic  part  of  the  SWT  Eq.(32)  is 
Af 

From  Rice  solution  Eq.(21)  and  Eq.(22)  one  obtains 


(40) 


As  _ 

*^max  2  ~ 


(l+n-) 

AK^ 

(l+n');TOoX 

^0 

M}.+Vl)kC^X 

d+n') 


(41) 


Substituting  RHS  of  Eq.(40)  into  Eq.(41)  and  after  rearrangement  the 
following  relationship  is  obtained 
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<TV«f'(2Nf)'=”  = 


-1 


4(l+n’) 

(l  +  n')nE^ 


1 

(1+n-) 


n2 


(42) 


Setting  X  =  p*  and  solving  for  Nf 


Nf=- 


-I 


40+n') 


(l  +  n');dEcr'f  '/9  * 


b+C  r- 


(K„ax)“"''(AK)<'«') 


2 

b+C 


(43) 


Where  Nf  is  a  number  of  cycles  required  to  break  the  process  zone 
element  and  to  advance  the  crack  by  the  ‘process  zone’  size  p*.  Therefore, 
the  corresponding  average  rate  of  crack  growth,  da/dN,  can  be  calculated 
as 


-1 

b+C 


da  _  p* 

dw"^ 


=  i2p*) 


40+n') 

(1  +  n');zEfT'j.  'p  * 


-2 

b+C 


(44) 


The  average  crack  growth  rate  in  the  Paris  region  can  be  expressed  in  the 
following  form 


(45) 


Where  AK*  — (AK^)^  ^(^max)^  parameter  driving  force 

proposed  by  Kujawski  [41]  and  Cp  and  mp  are  fitting  parameters. 
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Comparing  Eq.(44)  with  Eq.(45)  the  corresponding  exponents  are  related 
as 


n’ 

a  = - 

1  +  n' 


(46) 


A  similar  analysis  can  be  carried  out  for  the  near  threshold  region,  where 
the  crack  growth  is  dominated  by  elastic  stress  and  strain. 


Thus,  the  elastic  part  of  the  SWT  Eq.(32)  will  have  the  form  of 


The  stress  and  strain  at  the  process  zone  p*  can  be  calculated  using  the 
elastic  solution,  Eq.(31)  by  setting  x  =  p* 


Equating  RHS  of  Eq.(47)  and  Eq.(48)  and  solving  for  Nf  one  gets 


\(  K  AkV^ 
2  Anp  *  (jf  j 


(49) 


Considering  that  Nf  cycles  are  required  to  propagate  the  crack  by  the 
process  zone  size  p* 


The  average  crack  growth  rate  in  the  near  threshold  region  can  be 
expressed  in  the  following  form 
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(51) 


dN  J 

Comparing  Eq.(50)  with  Eq.(51)  and  Eq.(45),  the  following  relations 
result 


and 


(52) 


The  relations  (46)  and  (52)  demonstrate  that  the  exponents  of  fatigue  life 
and  crack  growth  equations  (or  the  slopes  of  the  corresponding  curves  in 
log-log  coordinates)  are  interrelated.  They  also  show  that  the  exponent  a 
in  the  two  parameter  driving  force,  AK*,  is  decreasing  from  Vi  at  the  near 


threshold  region  to  nV  (1+n’)  at  the  Paris  region.  This  indicates  that 
sensitivity  to  Kmax  is  decreasing  with  increasing  crack  growth  rate. 

In  the  following,  a  parametric  study  is  presented  that  illustrates  the 
sensitivity  of  the  relationship  between  the  slopes  of  the  SWT  curve,  SWT 
vs.  Nf,  and  the  crack  growth  rate  curve,  da/dN  vs.  AK,  in  both  elastic  or 
the  near  threshold  and  plastic  or  the  Paris  regions.  Fig. 25  shows  the 
relation  between  the  slopes  of  the  curves  in  the  elastic  region  of  the  SWT 
plot  versus  the  near  threshold  region  of  the  crack  growth.  The  relation 
corresponding  to  the  plastic  region  of  the  SWT  plot  versus  the  Paris 
region  of  the  crack  growth  is  depicted  in  Fig. 26.  Both  figures  indicate  the 
same  trend.  In  particular.  Fig. 25  indicates  that  the  variations  in  the  slope 
of  the  elastic  part  of  the  SWT  plot  affect  only  the  near  threshold  slope  of 
the  da/dN  vs.  AK  predictions.  Similarly,  Fig. 26  shows  that  for  the  slope 
variations  of  the  plastic  part  of  the  SWT,  the  corresponding  changes 
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occurred  only  in  the  Paris  region  slope  of  the  da/dN  vs.  AK.  It  can  be 
noted  that  somewhat  higher  sensitivity  is  observed  in  the  elastic  and  the 
near  threshold  regions  in  comparison  to  the  plastic  and  the  Paris  regions. 
This  is  consistent  with  the  Eq.(46)  and  Eq.(52)  for  the  Paris  and  the  near 
threshold  regions,  respectively. 

In  Figs. 25  and  26  some  fluctuations  in  the  predicted  curves  can  be  seen  at 
the  lower  part  of  the  Paris  region.  In  the  simulation,  two  straight  lines 
with  a  sharp  knee  point  were  used  instead  of  actual  smooth  SWT  curve. 
This  simplified  representation  of  SWT  curve  resulted  in  fluctuations  in 
the  predicted  curves  when  the  corresponding  lives  of  the  crack  tip 
elements  were  spread  around  the  sharp  knee  point. 
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rAE/2(1-R)  (Pa) 


7.  DETERMINATION  OF  THE  CRACK-TIP  ELEMENT  SIZE,  p* 

The  determination  of  the  crack  tip  element  size  or  the  process  zone  size, 
p*,  is  very  important  because  it  defines  the  distance  ahead  of  the  crack 
where  the  stress  profile  has  a  very  low  gradient  due  to  yielding  of  the 
material  and  crack-tip  blunting.  The  process  zone  size  has  been  taken  as 
constant  by  many  researchers  like  Glinka  [10-12],  Ellyin  and  Kujawski 
[15-18,  42]  and  Antolovich  [43].  Other  researchers  believe  that  the 
process  zone  size  is  variable  and  is  dependant  on  AK  [44]  or  it  may  be 
linked  to  the  microstructure  or  the  grain  size  of  the  material  [45]. 

In  this  model,  the  process  zone  is  considered  to  be  constant  for  a  given 
material.  The  process  zone  size,  p*,  is  calculated  using  one  point  of  the 
crack  growth  rate  data.  This  point  corresponds  to  the  transition  between 
the  Paris  and  near-threshold  regions  on  the  crack  growth  curve.  It  can  be 
taken  as  the  point  of  intersection  of  the  best  fit  lines  for  the  near¬ 
threshold  region  and  the  Paris  region.  Fig. 27  shows  a  simplistic  view,  as 
to  how  this  ‘knee-point’  is  determined  from  the  experimental  data. 

This  knee-point  should  satisfy  the  driving  force  equations  for  both  the 
Near-threshold  region  and  the  Paris  region.  By  comparing  these  two 
driving  forces  at  the  knee-point,  a  relation  for  p*  can  be  obtained. 
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Figure  27.  Best-fit  line  approach  to  determine  the  knee-point 
for  the  p*  determination. 


This  point  once  determined,  should  provide  the  same  da/dN  for  both  the 
equations  for  crack  growth  rate  in  Paris  region  and  the  near-threshold 
region. 

The  expression  for  crack  growth  rate  in  Paris-region  from  Eq.(44) 


-1 

-1 

b+c 

r 

da  p* 

4(l+n') 

n’  1 

(l  +  n');zEcr'f  'p  * 
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And  expression  for  da/dN  in  near-threshold  region  from  Eq.(50) 


f  n 


vNfy 


da 


=  2p- 


Atug 


fP*^ 


KLxAK2 


(54) 


Comparing  the  Eq.(53)  and  Eq.(54)  for  same  da/dN  at  the  knee-point 
following  equation  is  derived,  which  represents  the  final  expression  that 
is  used  to  calculate  p*. 
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8.  ACCUMULATION  OF  DAMAGE  AHEAD  OF  CRACK-TIP 
Damage  accumulation  in  the  elements  ahead  of  the  crack  tip  is  controlled 
by  the  Palmgren-Miner  linear  damage  rule.  Failure  life  for  each  element  is 
calculated  using  the  SWT  failure  criterion.  Damage  is  stored  in  the 
elements  within  the  plastic  zone  only,  as  it  was  seen  that  the  damage  past 
the  plastic  zone  is  negligible  and  can  be  considered  as  zero.  The  element 
is  considered  to  fail  when  the  damage  in  the  element  reaches  unity.  The 
detailed  explanation  of  damage  accumulation  process  and  the  element 
failure  is  presented  ahead  (see  section  9). 

The  important  thing  to  notice  is  the  gradient  of  the  damage  ahead  of  the 
crack  tip.  Fig. 28  and  Fig. 29  show  the  accumulated  damage  at  steady  state 
in  the  elements  ahead  of  the  crack  tip  for  various  Kmax  and  R-ratio  values. 
It  can  be  seen  in  Fig. 28  that  the  gradient  of  the  accumulated  damage  is 
very  high  close  to  the  crack  tip,  and  the  majority  of  the  damage  is  in  the 
elements  within  the  cyclic  plastic  zone.  After  the  cyclic  plastic  zone  the 
gradient  of  the  accumulated  damage  becomes  nearly  constant,  at  which 
point  the  damage  is  not  very  high. 
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Figure  28.  Accumulated  damage  distribution  in  elements  for 
4140  steel  (linear  scale). 


Another  important  thing  which  can  be  noted  from  Fig. 29  i.e.  the  log-log 
plot  of  damage  accumulation  is,  that  Kmax  controls  the  extent  of  the 
damage  zone  i.e.  the  monotonic  plastic  zone  size,  and  AK  values  control 
the  amount  of  damage  in  the  elements.  The  data  has  been  compared  for  R- 
ratio  0.1  and  0,  It  can  be  noted  that  for  the  same  K^ax  value,  accumulated 
damage  is  more  for  R  =  0. 
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Figure  29.  Accumulated  damage  distribution  in  elements  for  4140 
steel  (log-log  scale). 


As  explained  above,  Fig. 28  and  Fig. 29  show  the  fact  that  the  main  damage 
accumulation  ahead  of  the  crack  tip  is  in  cyclic  plastic  zone,  some  damage 
of  a  very  small  magnitude  is  there  up  to  the  monotonic  plastic  zone,  after 
that  the  damage  is  almost  negligible.  This  is  supported  by  the  statements 
in  references  [7-12,  15-18],  explained  on  the  basis  of  the  gradient  of  the 
stress-strain  profiles. 
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9.  CALCULATION  OF  THE  CRACK  GROWTH  RATE 
In  the  present  approach,  the  un-cracked  region  of  the  specimen  is  modeled 
as  a  set  of  uniaxial  specimens  of  width  equal  to  p*,  see  Fig. 24.  Then,  for 
a  given  load  cycle  the  maximum  and  minimum  stress-strain  distribution 
ahead  of  the  crack  tip  is  determined  using  the  Rice  and  elastic  solutions 
discussed  earlier. 

For  every  crack  tip  element,  the  SWT  parameter  term  .  is 

2(l-Rel) 

calculated 

Where  Rei  is  the  local  stress  ratio  of  amin/cTmax  for  that  particular  element. 
Based  on  the  product  for  the  SWT  parameter,  the  fatigue  life  Nf  for  each 
element  is  determined. 

Damage  stored  in  z*'’  element  is 


A 

Nf 


+  (previous  damage  ‘d”) 


(56) 


Where  n|-  =  Fatigue  life  of  element,  Nf=  Fatigue  life  for  the  z*'’ 
element. 


For  any  previously  damaged  element,  when  crack  tip  arrives  at  that 
element  the  remaining  fatigue  life  is  calculated  by  using  the  Palmgren- 
Miner  ‘linear  damage  rule’  [22],  and  by  keeping  track  of  how  much 
previous  damage  ‘d’  is  present  in  that  element. 

So  the  remaining  failure  life  for  an  element  is  calculated  as 

N,=(l-d)Nf  (57) 
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Basically,  this  gives  the  number  of  cycles  required  to  make  damage  equal 
to  one  in  that  element  at  crack  tip,  which  is  considered  broken  when  the 
damage  ‘d’  reaches  the  value  of  1. 

Similarly,  the  damage  for  the  elements  from  2"*^  element  onwards  can  be 
calculated  by  using  this  value  of  remaining  life  ‘N/  from  Eq.(57)  in  the 
following  way: 

.  _Nr  , 

°new-  j 
JN  f 

The  crack  growth  rate  da/dN  is  found  by  dividing  element  width  ‘p*’  by 
the  number  of  cycles,  Nr,  required  to  break  this  element  at  crack  tip. 

Given  by  equation: 


da  _  p* 


(59) 


This  process  is  repeated  for  every  element  and  the  crack  is  propagated  by 
one  element  each  time,  hence  accumulating  the  damage  in  the  elements 
ahead  of  the  crack  tip  and  breaking  the  elements  for  damage  equal  to 
unity. 
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Fig. 30  and  Fig. 31  ahead,  show  the  experimental  and  the  predicted  data 
using  the  above  discussed  damage  accumulation  approach,  Eq.(59),  for  the 
4140  steel  and  the  4340  steel  respectively. 


Figure  30.  Predicted  vs.  Experimental  crack  growth  data  for  4340 
steel,  at  R  =0.1  and  p*  =  1.56e'^. 
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Figure  31.  Predicted  vs.  Experimental  crack  growth  data  for  4140 
steel,  at  R  =0.1  and  p*  =  2.05e'^. 


It  can  be  noted  from  the  predictions  that  the  model  works  well  for 
constant  amplitude  loading.  The  p*  parameter  has  been  fitted  to  the 
experimental  data  for  this  prediction. 

Any  R-ratio  can  be  predicted  with  this  model,  condition  being  that  the 
experimental  data  should  be  known  before  hand,  so  that  p*  can  be 
determined. 

More  predictions  have  been  shown  for  various  materials  at  different  R- 
ratios  in  Appendix-A. 
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10.  OVERLOAD  EFFECTS  ON  CRACK  GROWTH  RATE 


The  accurate  prediction  of  fatigue  crack  growth  behavior  under  variable 
amplitude  loading  depends  on  accurate  modeling  of  load-interaction 
effects,  one  important  load  interaction  effect  being  the  retardation  of 
fatigue  crack  growth  due  to  application  of  an  overload  [48]. 

Extensive  investigations  have  been  conducted  [1-5  in  49]  in  an  attempt  to 
explain  the  phenomenon  of  crack  growth  retardation  resulting  from  tensile 
overloads.  These  investigations  have  resulted  in  proposal  of  several  crack 
growth  retardation  theories.  In  general,  these  theories  attribute  delayed 
retardation  behavior  to  crack  closure  [27],  residual  stresses  [49,  50,  51], 
or  a  combination  of  these  mechanisms  [52]. 

The  crack  closure  theory  suggests  that  the  delay  in  fatigue  crack  growth  is 
caused  by  the  formation  of  a  zone  of  residual  tensile  deformation  left  in 
the  wake  of  propagating  crack.  This  deformation  causes  the  crack  to 
remain  closed  during  a  portion  of  applied  tensile  load  cycle.  Hence,  the 
crack  extends  only  due  to  the  portion  of  the  load  cycle  that  is  above  the 
crack  opening  level. 

The  residual  theory  is  somewhat  similar  to  the  crack  closure  theory  except 
that  it  attempts  to  account  for  the  crack  growth  based  on  actual  forces 
acting  at  the  crack  tip.  It  focuses  more  on  what  is  happening  ahead  of  the 
crack  tip  and  at  the  crack  tip.  This  theory  suggests  that  the  application  of 
high  tensile  load  cycle  forms  residual  compressive  stresses  in  the  vicinity 
of  the  crack  tip,  and  these  stresses  reduce  the  rate  of  fatigue  crack  growth. 
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Lang  [51]  reports  that  the  major  restraint  to  crack  extension  is  caused  by 
residual  compressive  stresses  ahead  of  the  crack  tip.  Residual  compressive 
stresses  are  the  result  of  a  reverse  plastic  zone  created  after  a  remote 
tensile  stress  has  been  removed. 

In  this  study,  an  attempt  has  been  made  to  model  the  residual/internal 
stress  profiles  ahead  of  the  crack  tip  based  on  the  data  from  the  literature, 
to  explain  the  crack  growth  retardation  due  to  overload  and  effects  on 
crack  growth  due  to  block  loading. 
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10.1  MODELING  RESIDUAL/INTERNAL  STRESSES 


The  first  objective  for  this  kind  of  modeling  should  be  to  investigate  the 
distribution  of  residual  stresses  and  to  determine  the  correlation  of  these 
measurements  with  the  residual  stress  theories. 

Experimental  test  results  for  the  measurements  of  the  residual  stresses  at 
the  crack  tip,  before  and  after  application  of  the  overload  have  been 
analyzed;  from  the  studies  that  measured  the  stress  distributions  in  the 
vicinity  of  the  crack  tip  by  x-ray  diffraction  techniques  [49,  50  and  52]. 
Holloway  [49]  in  his  study  measured  the  residual  stress  profiles  before 
and  after  the  overload,  and  after  some  crack  extension  from  the  point  of 
application  of  the  overload  for  1020  Steel. 

The  residual  stress  profiles  before  and  after  the  overload  are  shown  in 
Fig. 32,  The  overload  is  twice  of  the  base  maximum  load. 

It  can  be  noted  from  the  Fig. 32  that  the  residual  stress  at  the  crack  tip  is 
quite  bigger  than  the  base  load,  the  point  of  minimum  stress  is  the  point 
of  application  of  overload. 
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-  -  -  PRIOR  TO  OVERLOAD 


Figure  32.  Residual  stress  profiles  measured  by  x-ray  diffraction  in 
1020  Steel  before  and  after  the  overload  [49]. 

Residual  stress  profile  after  overload  and  residual  stress  profile  after  a 
crack  extension  of  0.054  inches  past  the  overload  are  shown  in  Fig. 33. 
The  distribution  shows  that  the  residual  stress  profile  remains  essentially 
undisturbed  by  crack  extension. 
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POST  OVERLOAD 


Figure  33.  Residual  stress  profiles  after  overload  and  after  crack 
extension  of  0.054  inches  past  point  of  overload  [49]. 


Upon  some  more  crack  extension  to  0.111  inches  past  the  overload  point, 
the  maximum  residual  compressive  stress  was  found  to  be  reduced,  while 
the  stresses  behind  and  ahead  of  the  crack  tip  remained  essentially  the 
same,  Fig. 34  shows  the  measured  stress  profiles. 

Similar  trend  was  noted  when  the  crack  was  extended  further  to  0.265 
inches  past  the  overload  point.  The  point  of  overload  application  remained 
under  maximum  compressive  residual  stress,  while  there  were  some 
positive  stresses  seen  immediately  behind  the  crack  tip. 
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-  -  -  CRA®  TJP  8.0S4  INCH 


DISTANCE  AHEAD  OF  OVERLOAD  POINL  INCHES 

Figure  34.  Residual  stress  profiles  after  crack  extensions  of  0.054 
and  0.1 1 1  inches  past  the  point  of  overload  [49]. 


Results  from  this  suggest,  that  a  region  of  residual  compressive  stresses 
from  the  overload  exists,  and  is  evident.  This  region  of  residual 
compressive  stresses  remained  undisturbed  when  the  crack  advanced 
through  it. 

When  the  crack  tip  was  grown  past  the  influence  of  the  overload  induced 
plastic  zone,  tensile  stresses  were  found  in  the  wake  of  the  crack  between 
the  point  of  overload  application  and  crack  tip.  This  finding  is  in  direct 
contrast  with  the  crack  closure  concept  which  suggests  that  stresses 
behind  a  crack  tip  are  compressive  in  nature.  Holloway  [49]  concluded 
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that  overload  retardation  is  due  to  residual  compressive  stresses  generated 
due  to  the  overload  application,  and  that  the  compressive  stress  due  to 
this  overload  essentially  stays  undisturbed  by  the  crack  growth. 

Allison  [52]  conducted  a  similar  study  on  1045  steel  to  determine  the 
residual  stresses  before  and  after  the  application  of  overload. 

Fig. 35  shows  the  residual  stress  profiles  before  and  after  the  overload 
application  and  the  maximum  stress  profile  at  the  base  load. 


Figure  35.  Residual  stress  profiles  measured  by  Allison  [52],  before 
and  after  the  application  of  overload. 
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He  measured  the  residual  stress  profiles  after  a  crack  extension  of  1 .04 
mm  beyond  the  point  of  application  of  overload  and  is  shown  in  Fig. 36. 


Figure  36.  Residual  stress  profiles  after  the  overload  and  after  a 
crack  extension  of  1.04  mm  [52]. 

This  study  indicates  that  the  residual  stress  profile  behavior  before  and 
after  the  overload  and  after  the  crack  extension  is  very  similar  to  what 
was  shown  by  Holloway. 

Allison  also  reached  similar  conclusions  as  Holloway,  but  the  difference 
being  that  he  referred  to  closure  being  a  surface  phenomenon,  because  of 
the  stress  imbalance  that  exists  inside  and  outside  of  the  specimen. 
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The  analytical  modeling  of  the  residual  stress  can  be  done  with  the  Rice 
and  elastic  solution  at  the  unloaded  state  as  was  shown  earlier,  which  is 
consistent  with  the  model. 

Once  this  residual  stresses  have  been  modeled,  it  is  needed  to  apply  the 
effect  of  these  residual  stresses  on  the  driving  force.  For  this  Residual 
stress  intensity  factor  Kres,  or  Internal  stress  intensity  factor  Kint,  in  case 
of  block  loading  is  needed  to  be  calculated. 

To  calculate  the  Kres  or  Kint  from  the  stress  profiles,  weight  function 
method  has  been  used,  and  a  subroutine  has  been  developed.  The  detailed 
explanation  of  the  weight  function  procedure  and  its  adaptation  to  the 
current  model  is  given  in  the  next  section. 
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10.2  WEIGHT  FUNCTION  PROCEDURE 


The  Stress  Intensity  Factor  solutions  for  several  simple  cases  have  been 
published  in  the  handbooks,  but  for  several  non-uniform  stress  fields  and 
practical  cases,  the  handbook  solutions  are  inadequate. 

The  Weight  Function  method  proposed  by  Bueckner  [53]  and  Rice  have 
proved  to  be  very  useful  for  calculation  of  Stress  Intensity  Factors, 
especially  for  cracks  subjected  to  non-uniform  stress  fields.  Later  Glinka 
[54,  55]  and  Niu  [56,  57]  have  derived  weight  functions  for  various 
geometries  which  are  easy  to  be  incorporated  into  a  computer  algorithm 
for  quick  calculations. 

The  unique  feature  of  weight  function  approach  is  that  once  the  weight 
function  for  a  particular  cracked  body  is  determined,  the  stress  intensity 
factor,  K,  for  any  loading  system  applied  to  the  body  can  be  calculated  by 
a  simple  integration  of  the  expression  of  following  form: 

a 

K=  Jcr(x)m(x,a)dx  (60) 

0 

where  a(x)  is  the  stress  field,. m(x,a)  is  the  weight  function  for  the 
particular  geometry,  ‘x’  is  the  distance  from  the  crack  tip  and  ‘a’  is  the 
crack  length. 

The  technique  for  deriving  the  weight  function  was  proposed  by  Petroski 
and  Achenbach  [58],  and  later  by  Glinka  [54-57]. 

It  was  shown  by  several  authors  that  the  Petroski-Achenbach  method  can 
sometimes  lead  to  inaccurate  weight  functions.  The  crack  opening 
displacement  function  proposed  by  Petroski  and  Achenbach  for 
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determination  of  the  weight  functions  was  not  accurate  for  the 
discontinuous  and  high  gradient  reference  stress  fields. 


Glinka’s  weight  function  approach 

Ghnka  and  Niu  [54-57]  had  an  important  finding  that  the  weight  functions 
for  a  variety  of  crack  configurations  had  the  same  general  mathematical 
form.  They  derived  a  universal  weight  function  expression  which  could  be 
applied  to  a  wide  variety  of  Mode  I  cracks  regardless  of  geometry.  The 
associated  parameters  of  the  weight  function  corresponding  to  the 
particular  cracked  body  geometric  configuration  may  be  different. 

Glinka  et  al.  used  the  Petroski-Achenbach  approach  to  derive  a  general 
weight  function  of  the  form 


m(x,a)  = 


^27t{di-x) 


1  +  M, 


^  x^ 
1-- 


V 


-t-M- 


1-- 


(61) 


Where  m(x,a)  is  the  weight  function,  ‘a’  is  the  crack  length,  ‘x’  is  the 
distance  from  the  crack  tip,  and  Mi,  M2  are  the  weight  function 


parameters.  Details  about  the  determination  of  the  weight  function 
parameters  are  presented  in  the  discussion  ahead 

Accuracy  of  weight  function  approximated  by  Eq.(61)  depends  on  the  type 
of  the  reference  stress  intensity  factor  ‘K,.  and  the  associated  reference 


stress  field  ‘ar(x)’and  it  is  inaccurate  for  very  non-uniform  or 
discontinuous  local  reference  stress  field  (r,(x).  It  was  later  pointed  out  by 
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Glinka  and  Niu  that  in  some  cases  the  three-term  expression  of  Eq.(61) 
was  not  sufficient,  requiring  sometimes  more  than  three  terms. 


For  this  reason  a  general  weight  function  in  the  form  of  Eq.(62)  is  used 
2 


m(x,a)  = 


■yj2;r(a-x) 


1  +  Mi 


V  a. 


a; 


+  ....  +  M„  1-- 


(62a) 


m(x,a)  = 


^2K{a-x) 


1  +  M, 


1/2 


^  X^' 


.'-f) 


r 

1-- 
V  a; 


(62b) 


However,  it  is  not  known  how  many  terms  in  Eq.(62)  should  be  used  for  a 
particular  geometry  since  there  are  some  weight  functions  which  can  not 
be  approximated  by  using  even  seven  terms  in  the  expression. 

Although  other  forms  of  weight  functions  are  admissible  for  crack 


problems,  expressions  (61)  and  (62)  are  the  most  frequently  used. 

To  determine  the  stress  intensity  factor,  weight  function  method  requires 
an  integration  of  the  stress  field  and  the  weight  function,  where  the 
weight  function  can  be  evaluated  from  the  expression  for  Universal 
weight  function  of  the  form  of  Eq.(61)  and  Eq.(62). 


^^•2-2  Determination  of  the  weight  function  for  a  eeometrv 

The  technique  that  has  to  be  followed  to  evaluate  the  constants  Mi,  M2, 

M3  ,...Mn  is  explained  below  which  has  been  proposed  by  Glinka. 

It  was  shown  by  Glinka  [54,  55]  that  the  three  term  expression  of  Eq.(61) 
can  approximate  weight  functions  for  shallow  edge  and  central  through 
cracks  of  relative  depth  a/t  <  0.5,  while  the  four-term  form  of  Eq.(62)  can 
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be  sufficiently  accurate  for  a  variety  of  cracks  without  any  limitations 
regarding  their  depth. 


10-2.2.1  General  approach  for  ID  cracks 


If  the  four  term  expression  of  Eq.(62)  is  used,  the  derivation  of  a  weight 
function  for  any  particular  geometry  can  be  reduced  the  determination  of 
three  parameters  Mi,  M2,  and  M3.  Thus,  the  weight  function  m(x,a)  can  be 
determined  without  crack  opening  displacement  function.  Only  two 
reference  stress  intensity  factors  have  to  be  known.  Once  these  are  known, 
a  set  of  three  equations  can  be  defined  for  the  three  reference  stress 
intensity  factors  Kri,  Kr2  and  Kr3. These  three  equations  involve  the  three 
constants  Mi,  M2,  and  M3.  So  the  equations  can  be  derived  by  integration 
ofEq.(60) 


If  the  four  term  expression  of  Eq.(62b)  is  used. 


1^  =  'f  ^^ri(x) 


1  +  Mi  1- 


a; 


-fM- 


1-^ 

V  a; 


+  M, 


1-^ 
k  a; 


dx  (63a) 


K  =  "r  20-^2 (x) 
Q^j2n(a-x) 


1  +  M, 


1-- 

V  a; 


+  M- 


^  x^ 
1-- 

k  a; 


■fM. 


1-^P 


dx  (63b) 


-  'f  ^<^r3(x) 


1  +  M,  1-- 
k  a; 


+  M. 


x^ 

1-- 

a; 


+  M. 


1- 

V  a 


X  \2 


dx  (63c) 
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Eq.(63a)  -  Eq.(63c)  can  be  used  as  a  set  of  three  simultaneous  algebraic 
equations  which  can  be  solved  for  Ml,  M2  and  M3  and  these  values  can 
be  used  to  calculate  SIF  associated  with  any  stress  field. 


10.2.2.2  General  approach  for  2D  cracks 

For  two-dimensional  cracks  it  is  necessary  to  calculate  the  SIF  for  at  least 
two  points 

The  deepest  point,  A 

a 

Ka=  Jo-A(x)mA(x,a)dx  (64) 

0 

The  surface  point,  B 

a 

Kb=  /(TB(x)mB(x,a)dx  (65) 

0 

The  weight  functions,  mA(x,a)  and, mB(x,a)  are  calculated  as  follows: 


mA(x,a)  = 


1  +  M 


lA 


1-^ 

a; 


+  M 


2a: 


l-2i 
.  ay 


+  M 


3AI 


1-- 
.  a 


X  12 


(66) 


mB(x,a)  = 


1  +  M 


IB 


1  X  2 

1--  +M 


2B 


^  X^ 
1-- 

V  a; 


+  M 


(67) 


The  parameters  MIa,  M2a,  M3a  for  the  deepest  point,  and  MIb,  M2b,  M3b 
for  the  surface  point  are  calculated  depending  on  the  crack  geometry. 
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^0.2.3  Numerical  algorithm  for  weight  function  procedure 

The  SIF  solutions  by  weight  function  method  require  huge  functions  to  be 

integrated,  so  analytical  integration  requires  too  much  calculation  time 

and  can  be  very  rarely  solved  due  to  its  complexity.  So  an  algorithm  based 

on  Glinka’s  weight  function  has  been  followed  to  solve  the  integration 

numerically. 


The  general  weight  function  expression  in  the  form  of  Eq.(60)  can  be 
written  as 


K=  J- 


cr(u)2a 


1  +  M,(1-u)2  +M2(1-u)  +  M3(1-u): 


du 


(68) 


Where  u  =  x/a 

Then  the  integration  domain  0  <u  <1,  may  be  divided  into  “n”  sub¬ 
intervals,  in  each  of  which  the  stress  function  a  (u)  is  approximated  by  a 
straight  line  similar  to  treatment  shown  in  Fig. 37. 


cri(u)  =  AjU  +  Bi  ,  Ui_i<u<Ui 


(69) 


where 


_<^(Ui)-g'(Ui-,) 


Ui-Ui_i 


Bi  =  (T(u,)-AiUi 


(70) 


The  number  of  sub  intervals  must  be  chosen  with  such  a  way  as  to 
approximate  the  stress  function  by  straight  lines  with  reasonable  accuracy. 
So  now  the  Eq.(68)  can  be  written  in  the  form  of  Eq.(71)  below 


-i-i 

K  =  a  Jm(u,a)(AiU  +  Biu)lu 


(71) 
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where 


m(u,a)  =  1  +  M,(l  -  u)i  +  M2(1  -  u)+ M3(1  -  u)| 

So  finally  the  stress  intensity  factor  can  be  calculated  from 

K  =  a^S*m(a,U*] 
i=l 


(72) 


(73) 
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Figure  37.  Illustration  of  integration  method  based  on  piecewise 
linearization  of  the  stress  function  a)  Weight  function  m(x,  a), 
b)  Piecewise  stress  function. 
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The  term  m(a,U*i )  represents  the  value  of  the  weight  function  for  u  =  U*i 
The  parameter  S  j  and  U  i  are  the  area  under  the  normalized  stress  curve 
<7  (u)  and  the  coordinate  ‘u’  of  its  centroid,  respectively,  corresponding  to 
the  sub  interval  “i”  similar  to  that  shown  in  Fig. 37. 

S- =^(CT(Uj)-o-(Uj_,)Xui-Ui_,)  (74) 


[2<T(Ui_i)  +  cr(Ui)] 

3[cr(Ui)  +  a(Ui_,)] 


(75) 


The  singularity  of  the  weight  function  near  the  end  of  the  integration 
interval  “i”  where  (u^  1)  makes  it  inappropriate  to  use  Eq.(73)  in  this 
region,  for  0.95  <u  <1.  Therefore  for  this  particular  interval  or  for 


some  very  other  small  interval  depending  upon  the  stress  profile 
analytical  expression  can  be  derived  for  the  integration.  The  analytical 
expression  can  be  used  for  the  singularity  condition  in  the  algorithm  that 
has  been  used  for  the  software,  summation  of  the  analytical  expression  is 
used  throughout  the  whole  crack  length 


The  results  for  different  geometries  were  matched,  for  the  SIF  calculation 
by  Glinka’s  FALPR05  application  and  the  developed  MATLAB  procedure. 
The  results  matched  very  accurately  with  differences  of  less  than  about 
0.1%.  The  calculation  time  required  has  been  reduced  using  analytical 
piecewise  integration  expression. 

The  comparison  graphs  of  some  results  have  been  shown  Fig. 38  -  Fig. 41: 
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S.I.F  (MPaVm) 


-500.0 


Crack  length  (mm) 


Figure  38.  Typical  residual  stress  profiles  ahead  of  the  crack  tip 


Crack  length  (mm) 


Figure  39.  Comparison  of  S.I.F  calculation  for  the  stress  profiles 
by  Matlab  procedure  and  Glinka’s  FALPR05 


There  is  around  0.1%  difference  in  the  results  calculated,  which  arises  due 
to  the  analytical  piecewise  integration  used  throughout  in  the  Matlab 
procedure  to  boost  the  calculation  speed. 

To  see  how  the  good  the  procedure  works  for  very  sharp  gradient  profiles, 
Residual  stress  profiles  in  Fig. 40  are  the  data  digitized  for  situation  after 
a  spike  overload  at  the  interior  and  at  the  surface  of  the  specimen. 


Figure  40.  Residual  stresses  ahead  of  the  crack  tip  after  application 
of  a  spike  overload  at  interior  and  surface  of  the  centre  crack 
specimen. 
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Fig.41  below  shows  the  calculated  S.I.F  for  the  stress  profiles  after  the 
spike  overload  using  FALPR05  and  Matlab  procedure. 


Figure  41.  S.I.F  calculation  by  FALPR05  and  Matlab  procedure  for 
the  residual  stress  profiles  in  a  centre  crack  specimen  after  a 
spike  overload. 


The  difference  in  the  calculation  between  Glinka’s  FALPR05  and  the 
Matlab  procedure  is  very  small,  as  it  can  be  seen  from  the  comparison 
study  performed  above.  For  even  very  sharp  gradient  profiles,  the 
accuracy  of  the  procedure  does  not  vary.  So  this  procedure  can  be  utilized 
effectively  to  evaluate  the  S.I.F  by  weight  function  procedure. 
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10.3  RESIDUAL  STRESS  PROFILE  MODELING 

For  the  residual  stress  prediction  ahead  of  the  crack  tip,  rice  and  elastic 
solution  in  unloading  have  been  used.  The  minimum  stress  profile  at  the 
zero-load  is  considered  as  the  residual  stress  in  the  material. 

The  results  from  Holloway’s  [49]  x-ray  diffraction  measurement  study 
were  modeled  after  unloading  to  zero  from  base  load  and  overload.  Fig. 42 
shows  the  prediction  and  the  experimental  data.  The  material  used  is  1020 
steel,  the  base  load  maximum  stress  intensity  factor  is  25.38  MPaVm  and 
the  overload  maximum  stress  intensity  factor  is  50.774  MPaVm. 


Figure  42.  Residual  stresses,  experimental  vs.  modeled,  before 
and  .after  the  overload  application. 
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There  is  some  similarity  in  the  trends  in  Fig. 42,  but  close  to  the  crack  tip 
the  stresses  observed  by  the  experimental  procedure  are  very  less, 
compared  to  those  predicted  using  Rice  and  elastic  solution.  This  can  be 
attributed  to  the  fact  that  the  x-ray  diffraction  method  uses  the  average 
value  of  the  stress  within  the  region  at  which  the  beam  is  focused.  The 
size  of  the  x-ray  beam  spot  used  in  experimental  study  by  Holloway  was 
0.635  mm  in  diameter  on  the  specimen  surface. 

There  are  various  approaches  that  have  been  tried  to  model  the  residual 
stress  intensity  factor,  K,es,  profile  that  will  be  affecting  the  applied 
stress  intensity  factor  in  the  overload  affected  region.  The  logical  ones 
are  (a)  to  use  the  residual  stress  profile  after  the  overload  to  calculate  the 
residual-K  profile  or  (b)  to  use  the  difference  between  the  minimum  stress 
at  the  base  load  and  the  residual  stress  after  the  overload  to  calculate  the 
residual-K  profile. 
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10.4  PREDICTIONS  WITH  THE  MODEL 

The  retardation  of  the  overload  is  predicted  for  the  experimental  data 
from  Holtz  [46],  and  the  analysis  is  presented  ahead. 

The  data  has  been  obtained  for  4140  Steel  at  AK  =  19.3  MPaVm  at  R  =0.1, 
and  the  K^ax  overload  is  twice  that  at  the  base  load.  Fig. 43  shows  the 
stress  distribution  before  and  after  the  overload. 


Figure  43.  Stress  profiles  ahead  of  the  crack  tip  at  base  load  and 
after  overload  for  4140  steel. 


Residual  stress  after  the  overload  is  considered  to  stay  in  the  material 
ahead  of  the  crack.  Based  on  this  residual  stress  profile,  Kres  is  calculated, 
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which  will  affect  the  applied  K„ax  and  the  applied  AK.  Fig. 44  shows  the 
Kres  profile  for  the  modeling. 


Figure  44.  K  residual  profile  calculated  from  the  modeled  residual 
stress  in  the  material,  after  the  overload  application  by  using  weight 
function  approach. 

The  residual  K  affects  the  applied  Kmax  and  AK  in  the  following  way; 


^max 

(76) 

forK„i„>0 

(77) 

K^min  = 

for  K^i„<0 

(78) 

96 


Fig. 45  shows  the  predicted  crack  growth  rate  obtained  by  this  approach 


Figure  45.  Effect  of  overload  modeled  for  4140  steel  data  from 
Holtz  [46],  using  the  Kres. calculated  by  weight  function  method. 


From  Fig. 45  it  can  be  seen  that  the  retardation  obtained  by  the  model  is 
more  than  the  retardation  obtained  by  experiments.  The  maximum 
retardation  point,  in  terms  of  crack  length,  obtained  by  the  model  is  very 
close  to  the  one  obtained  by  experiments.  The  residual  K  profile 
calculation  needs  some  refinement  to  predict  the  real  behavior. 

In  the  second  approach,  residual  K  is  calculated  from  the  stress  profile 
obtained  from  the  difference  between  the  residual  stress  at  overload  and 
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the  minimum  stress  at  base  load,  assuming  that  the  overload  ‘over-writes 
the  base  load  minimum  profile.  Fig. 46  shows  the  stress  profile. 


Figure  46.  Residual  stress  profile  obtained  from  difference  between 
the  residual  stress  at  overload  and  minimum  stress  at  base  load. 


The  K  residual  profile  has  been  calculated  from  the  stress  profile  in 
Fig. 46  using  the  weight  function.  The  predicted  crack  growth  data  for  the 
overload  affected  region  is  shown  in  Fig. 47.  It  can  be  seen  from  this 
prediction  that  there  is  a  longer  retardation  in  the  predicted  crack  growth 
rate.  This  approach  of  using  the  difference  between  the  residual  stress  at 
overload  and  minimum  stress  at  base  load,  based  on  the  assumption  that 
the  minimum  stress  profile  at  base  load  stays  in  the  material  as  a  steady 
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state  affecting  the  driving  force.  Fig. 47  shows  that  the  method  explained 
above  does  not  give  good  results  in  comparison  to  the  approach  used  for 
obtaining  results  in  Fig. 45. 


Figure  47.  Crack  growth  retardation  due  to  overload  predicted  using 
the  residual  K  approach  for  4140  Steel. 


Another  set  of  overload  data  from  Stoychev  [59]  is  modeled  using  the 
approach.  For  this  data  the  base  load  was  applied  at  R=0.1  and  a  single 
overload,  twice  the  base  load  was  applied  and  the  crack  growth  rate  was 
acquired  using  the  optical  measurements.  Fig. 48  shows  the  3  sets  of 

experimental  data  and  the  predicted  data  using  K^es  for  2324  Aluminum 
alloy. 
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Figure  48.  Overload  data  for  2324  Al.  alloy  from  [59],  predicted 
using  Kres  approach. 


The  crack  growth  retardation  predicted  by  the  model  is  more  than  the 
experimental  data.  The  above  prediction  is  based  on  the  constants  found 
by  fitting  of  the  experimental  data  for  the  constant  R=0.1  for  2324-T39 
alloy.  Fig. 49  shows  the  fitted  constant  R-ratio  data.  The  slope  used  in  the 
near  threshold  region  was  mth=15.6.  Since  there  is  not  much  experimental 
data  in  the  near  threshold  region  the  best  fit  slope  was  estimated. 

In  the  overload  affected  region  the  AK  is  close  to  the  values  in  the  near 
threshold  region.  So  the  prediction  for  the  crack  growth  retardation 
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overload  should  be  very  sensitive  to  the  slope  of  crack  growth  curve  in 
the  near  threshold  region. 


Figure  49.  Constant  R-ratio  data  fitted  to  find  the  knee  point  and 
the  constants  ‘b’  and  ‘c’  for  2324-T39  Al.  alloy  at  R=0.1. 


To  verify  the  observation,  that  the  more  retardation  predicted  by  the 
model  may  be  due  to  the  high  slope  used  in  the  near  threshold  region  the 
experimental  data  was  again  fitted  with  a  lesser  slope  in  the  threshold 
region  with  the  knee  point  being  the  same.  New  values  for  constants  ‘b’ 
and  ‘c’  were  determined  and  the  simulation  was  run  again  for  the  crack 
growth  retardation  due  to  single  overload  application. 
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Fig. 50  shows  ihe  fitted  constant  load  ratio  data  with  a  slope  mth=l  1 .4  in 
the  near  threshold  region.  Based  on  the  constants  determined  from  this 
fitting  the  results  for  the  prediction  are  shown  in  the  Fig.51. 


Figure  50.  Constant  R-ratio  data  fitted  with  a  slope  of  11.4  in  the 
near  threshold  region  to  see  the  sensitivity  on  the  prediction. 
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Figure  51.  Sensitivity  of  the  model  to  the  slope  in  the  threshold 


region  to  predict  the  crack  growth  rate  in  overload  affected  region. 

It  can  be  noted  that  for  a  change  in  the  slope  mth  from  15.6  to  11.4,  the 
minimum  crack  growth  rate  for  the  overload  data  shifted  by  order  of  more 
than  1  decade  on  the  crack  growth  rate  axis.  The  model  is  very  sensitive 
to  the  slope  in  the  near  threshold  region.  So  for  better  prediction  results 
with  the  model,  there  should  be  sufficient  experimental  data  in  the  near 
threshold  region.  ' 
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11,  PREDICTION  OF  THE  R-RATIO  EFFECT  WITH  THE  MODEL 
Load  ratio  effect  was  predicted  using  the  current  model.  A  constant  value 
of  p*  was  used  to  predict  the  R=0.1  and  R=0.5  for  4340  steel.  The 

experimental  data  is  taken  from  [15].  Results  for  the  prediction  are  shown 
in  the  Fig. 52. 


Figure  52.  Experimental  and  predicted  crack  growth  rate  data  for 
4340  steel,  R  =0.1  and  0.5. 


The  developed  model  does  not  predict  the  R-ratio  effect  very  well.  There 
is  a  very  small  shift  in  the  crack  growth  rate  for  different  R-ratios.  This 
can  be  related  to  the  fact  that  this  model  uses  the  local  load-ratio  in  the 
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elements  to  predict  for  the  failure  life  from  the  SWT  equation.  Since  very 
close  to  the  crack-tip  the  local  R  «-l  or  has  a  very  low  gradient  and  is 
nearly  constant.  The  local  R-ratio  in  the  elements  ahead  of  the  crack  tip 
modeled  by  rice  and  elastic  solution  is  shown  in  Fig. 53. 


Figure  53.  Local  R-ratio  in  the  elements  ahead  of  the  crack  tip. 


Since  the  maximum  damage  is  present  in  the  elements  within  the  cyclic 
plastic  zone.  The  local  R-ratio  of  these  elements  controls  the  predicted  R- 
ratio  shift  by  the  SWT  failure  criterion. 

To  predict  the  R-ratio  effects  with  the  current  model,  p*  can  be  used  as  a 
fitting  parameter  to  fit  the  crack  growth  rate  to  the  R-ratio  to  be  predicted. 
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The  data  for  Fig. 52  is  fitted  with  a  different  p*  for  each  R-ratio,  the 
predicted  results  are  shown  in  Fig. 54. 


Figure  54.  Crack  growth  rate  predicted  for  R=0.1  and  0.5  for  4340 
steel  using  different  p*  value  for  each  R-ratio. 


The  different  p*  to  fit  the  data  for  each  R-ratio  works  well,  but  the 
experimental  data  should  be  known  before  hand.  So  the  proper  modeling 
of  the  local  R-ratio  in  the  first  few  elements  is  very  important. 


12.  CONCLUSIONS 


As  mentioned  earlier,  the  main  objective  of  this  study  was  to  develop  a 
mathematical  model  for  the  fatigue  crack  propagation  utilizing  the  stress 
and  strain  state  near  the  crack  tip  region  and  the  damage  accumulation 
process.  The  model  has  been  developed  and  a  Matlab  based  application  is 
made  to  run  all  the  simulations  and  obtain  the  results.  The  study  leads  to 
the  following  conclusions: 

(a)  The  Rice  solution  for  the  local  stresses  and  strains  ahead  of  the 
crack  tip  within  the  plastic  zone  and  elastic  solution  for  region 
outside  the  plastic  zone  works  good  with  the  model,  and  to  make 
the  distribution  continuous,  Irwin  correction  for  the  crack  tip  needs 
to  be  applied  to  the  elastic  solution. 

(b)  The  fatigue  exponents  in  the  SWT  equation  have  been  related  to  the 
slope  of  the  crack  growth  curve  in  the  Paris  and  near  threshold 
region. 

(c)  The  approximate  value  of  the  fatigue  exponents  ‘b’  and  ‘c’  for  a 
material  can  be  determined  from  the  average  of  the  slopes  of  the 
crack  growth  rate  curves  for  different  R-ratios. 

(d)  Both  Kmax  and  AK  are  the  contributing  parameters  to  the  driving 
force  for  the  model  which  is  consistent  with  the  recently  proposed 
two  parameter  approach. 

(e)  Constant  R-ratio  data  for  various  materials  has  been  predicted  with 
the  current  model,  with  process  zone,  p*,  being  a  calibrating 
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parameter.  The  predictions  are  in  good  agreement  with  the 
experimental  data. 

(f)  The  damage  accumulation  with  local  R-ratio  underestimates  the 
load  interaction  effects  on  the  crack  growth  rate,  Variation  in  the 
damage  and  the  local  R-ratio  does  not  predict  enough  retardation  or 
acceleration  as  predicted  by  the  experimental  data. 

(g)  Residual  stresses  are  assumed  to  stay  in  the  material  due  to  the 
overload  application;  these  residual  stresses  ahead  of  the  crack 
have  been  used  in  addition  to  the  variation  in  the  local  R-ratio  and 
accumulated  damage  to  model  the  overload  effects.  Weight 
function  methodology  has  been  used  to  calculate  Kres- 

(h)  The  results  of  the  Kres  calculation  from  the  four  term  weight 
function  expression  have  been  compared  to  the  results  from  the 
FALPR05  application  by  Glinka  and  the  results  match  very  well. 

(i)  It  has  also  been  shown  in  the  study  that  the  amount  of  retardation 
is  very  sensitive  to  the  slope  of  the  elastic  part  of  the  SWT  curve 
or  the  slope  of  the  crack  growth  curve  in  the  near  threshold  region. 

(j)  The  predicted  results  for  the  crack  growth  retardation  are  in 
reasonable  agreement  to  the  experimental  data  and  are  considered 
qualitatively  valid  for  overload  effect  in  different  materials. 
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13.  RECOMMENDATIONS  FOR  FUTURE  WORK 
Working  on  this  model  based  on  damage  accumulation,  some  limitations 
in  its  application  were  found  and  explored  while  many  other  questions 
were  developed  as  a  result. 

(a)  It  is  believed  that  through  some  future  exploration  p*  parameter 
can  be  made  independent  of  the  load  ratio  to  be  a  material 

parameter  rather  than  being  a  calibrating  parameter  for  each  load 
ratio. 

(b)  Study  should  be  done  to  improve  the  local  R-ratio  prediction  within 
the  cyclic  plastic  zone  with  the  model.  It  is  believed  that  local  R- 
ratio  within  cyclic  plastic  zone  is' not  very  realistic  which  might  be 
a  reason  for  the  inability  of  the  model  to  predict  the  R-ratio  effect. 

(c)  There  is  room  for  improvement  in  the  prediction  of  the  load 
interaction  effects;  future  work  should  be  focused  on  the  numerical 
determination  of  the  stress-strain  distribution  and  better  estimation 
of  the  parameters  influencing  the  crack  tip  driving  force. 
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APPENDIX  A 


KNEE-POINT  DETERMINATION  TO  FIND  p*  AND 
CONSTANT  R-RATIO  PREDICTIONS 
FOR  DIFFERENT  MATERIALS 
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da/dN  (m/cycle) 


Figure  Al.  Knee  point  used  to  find  p*  for  R  =0.1,  2324-T39  Al.  alloy  [59]. 


Figure  A2.  Predicted  vs.  experimental  data  for  R=0.1,  2324-T39  Al.  alloy  [59]. 
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Figure  A5.  Knee  point  used  to  find  p*  for  R  =0.1,  7475  Al.  alloy  [59]. 


Figure  A6.  Predicted  vs.  experimental  data  for  R=0.1,  7475  Al.  alloy  [59]. 


118 


da/dN  (m/cycle) 
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AK  (MPaVm) 

Figure  A7.  Knee  point  used  to  find  p*  for  R  =0.1,  A533-B1  steel  [15]. 


Figure  A8.  Predicted  vs.  experimental  data  for  R=0.1,  A533-B1  steel  [15], 
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Figure  A9.  Rnee  point  used  to  find  p*  for  R  =0.3,  A533-B1  steel  [15], 


Figure  A 10.  Predicted  vs.  experimental  data  for  R=0.3,  A533-B1  steel  [15]. 
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da/dN  (m/cycle)  ^  (m/cycle) 


da/dN  (m/cycle) 


Figure  A13.  BCnee  point  used  to  find  p*  for  R  =0,  0050A  cast  steel  [15], 


Figure  A14.  Predicted  vs.  experimental  data  for  R=0,  0050A  cast  steel  [15], 
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Figure  A17.  Knee  point  used  to  find  p*  for  R  =0,  C-Mn  cast  steel  [15]. 
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Figure  A18.  Predicted  vs.  experimental  data  for  R=0,  C-Mn  cast  steel  [15]. 
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Figure  A19.  Knee  point  used  to  find  p*  for  R  =0.5,  C-Mn  cast  steel  [15]. 


AK  (MPaVin) 


Figure  A20.  Predicted  vs.  experimental  data  for  R=0.5,  C-Mn  cast  steel  [15]. 


da/dN  (m/cycle) 


Figure  A21.  Knee  point  used  to  find  p*  for  R  -0,  Mn-Mo  cast  steel  [15], 


Figure  A22.  Predicted  vs.  experimental  data  for  R=0,  Mn-Mo  cast  steel  [15], 
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da/dN  (m/cycle)  asJ6N  (m/cycle) 
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'igure  A23.  Knee  point  used  to  find  p*  for  R  =0.5,  Mn-Mo  cast  steel  [15]. 


Figure  A24.  Predicted  vs.  experimental  data  for  R=0.5,  Mn-Mo  cast  steel  [15]. 
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da/dN  (m/cycle) 


Figure  A25.  Knee  point  used  to  find  p*  for  R  =0,  8630  cast  steel  [15]. 


Figure  A26.  Predicted  vs.  experimental  data  for  R=0,  8630  cast  steel  [15]. 
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igure  A27.  Knee  point  used  to  find  p*  for  R  =0.5,  8630  cast  steel  [15]. 
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Figure  A28.  Predicted  vs.  experimental  data  for  R=0.5,  8630  cast  steel  [15] 


APPENDIX  B 

WEIGHT  FUNCTION  CONSTANTS  FOR  DIFFERENT  GEOMETRIES 
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The  crack  geometries  and  the  associated  expressions  for  Mi,  M2  and  M3 
used  for  the  weight  function  calculation  are  given. 


Center  crack  specimen 

Mj  =  0.06987  +  0.401 17a-  5.5407a^  +  50.0886a^  -  200.699a'‘ 

+  395.5520^ -377.939a^  +  140.218a’ 

M2  =  -0.09049  -  2.14866a+  22.53250^  -  89.6553qP  +  210.599a'' 
-239.4450^ +11 1.1 28a^ 

M3  =  0.24721  + 2.56001a- 29.6349a’  + 1 38. 4a^  -347.255a'' 

+  457.128a^  -295.882a^  +  68.1575a’ 

where  a  =  —  =  crack  length 
w  specimen  width 


Double-edge  crack  specimen 

Ml  =0.08502 -0.02230a- 1.41028a’ +4.645590^  +  19.6924a'' 
-148.2660^ +336.837a^ -336.591a’ +127.009a^ 

M2  =  0.2234  -  0.61 46a+ 1 1 . 1 687a’  -  56.5326a’  +151 .937a'' 
-182.6340^+ 86.473  la^ 

M3  =  0.4983  +  0.75121a-10.5597a’  +  47.9251a’  -1 15. 933a'‘ 

+  131.9760^ -59.8893a^ 
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where 


a  = 


a 

w 


crack  length 
specimen  width 


Edge  crack  specimen 

(-  0.029207  +  a(0.2 1 3074  +  a(-  3.02955  +  a(5.901 933  +  a(-  2.6578))))) 

'  (l  +  a(- 1 .2597  +  c{-  0.048475  +  a(0.48125  +  c{-  0.5267  +  a(0. 34501)))))) 

jyf  (0.451 16  + 0(3.4624 +  a(-l  .07846 +  q(3.55857  +  q(- 7,55353))))) 

^  (l  +  o(- 1 .4966  +  a(0.7646  +  o(-  0.6593 1  +  a(0.2585  +  a(0. 1 14568)))))) 

M  -  (0-427195 +  a(- 3.7301 14 +  o(l6.2763  +  o(- 18.7999  +  0(14. 11211))))) 

^  (l  +  «(- 1 . 12919  +  o(0.033758  +  o(0.1921 14  +  o(-  0.65824  +  a(0.55466)))))) 

where  a  =  —  =  crack  length 
w  specimen  width 
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APPENDIX  C 
MATLAB  SUBROUTINES 
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Subroutine  to  calculate  the  p*  parameter 


% - Executes  on  button  press  in  ca.lc_dstar. 

function  calc_dstar_Callback (hObject ,  eventdata,  handles) 

%  hObject  handle  to  calc_dstar  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 


E  =  str2double (get (handles. editl, 'String' )) *1000000000; 
sy  =  str2double (get (handles. edit2, 'String' )) *1000000; 
n  =  str2double (get (handles. edits, 'String') ) ; 

SF  =  str2double (get (handles . edit4 , 'String' )) *1000000 ; 
eps  =  str2double (get (handles . edits, ' String ')) ; 
c  =  str2double (get (handles . edits, ' String' )) ; 
b  =  str2double (get(handles. edit?, 'String' )) ; 

Die  =  str2double (get (handles. editl4, 'String' )) *1000000; 
dadnch  =  str2double (get (handles .editlS, 'String')); 


setappdata (gebf , 
setappdata (gebf , 
setappdata (gebf, 
setappdata (gebf, 
setappdata (gebf , 
setappdata (gebf, 
setappdata (gebf , 


'E' ,  E) ; 
'sy,  sy)  ; 
'n',  n) ; 
'SF',  SF) ; 
'eps',  eps) ; 
'C,  c)  ; 

'b',  b) ; 


%  Units  = 
%  Units  = 

%  Units  = 


%  Units  = 


Pa 

Pa 

Pa 


Pa 


%  getting  applied  R  Ratio  from  applied  condition 

Ratio  =  str2double (get (handles .editl3 , 'String' )) ; 
kmax  =  (Dk/ (1-Ratio) ) ; 
kmin  =  Dk-kmax; 
clear  Ratio; 


%  MY  METHOD  FOR  FINDING  DSTAR  BY  USING  THE  KNEE  POINT  RELATION  THAT 
%  BOTH  ELASTIC  SWT  AND  PLASTIC  SWT  ARE  SATISFIED  THERE,  SO  EQUATING 
%  BOTH  EQUATIONS  FOR  CONSTANT  da/dN  RHO*  CAN  BE  FOUND  OUT. 


expl  ■=  (l/(i+n)),v 

exp2  =  (n/(l+n)); 
exp3  =  (l/(b+c)); 
exp4  =  (2/(b+c)); 
exp5  =  (1/ (2*b) ) ; 
exp6  =  (1/ (exp5-exp3) ) ; 

prtl  =  ( ( (4^expl) * (1+n) *pi*E*SF*eps) *exp3) ; 
prt2  =  ( ( (kmax*Dk) / (4*pi* (SF^2) ) ) *exp5) ; 
prt3  =  ( ( (kmax*exp2) * (Dk^expl) ) *exp4) ; 

Dstar  = ( ( (prtl*prt2) / (prt3) ) "exp6) ; 

% - 


setappdata (gebf ,  'Dstar ', Dstar) ; 
dstar  =  num2str (Dstar) ; 
set (handles. editie, ' String ', dstar) ; 
set (handles. editie, 'enable' , 'on' )  ; 
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Subroutine  for  the  calculation  and  plotting  of  the  SWT  table 


%  -  Executes  on  button  press  in  SWT_calc. 

function  SWT_calc_Callback (hObject ,  eventdata,  handles) 

%  hObject  handle  to  SWT_calc  {see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  GENERATING  SWT  TABLE 

%  In  the  SM  table  1st  Column  should  be  the  Product- (min.  to  max.),  2nd 
%  column  should  be  Nf  (max.  to  min.) . . . .Input  from  user  can  be  taken  in 
any  way 


%  MAKING  SM  TABLE  BY  USING  KNEE  POINT  AND  ELASTIC  AND  PLASTIC  PART 

Ds'tar  =  getappdata  (gcbf ,  'Dstar  * )  ; 
n  =  getappdata (gcbf , 'n' ) ; 
sy  =  getappdata (gcbf sy' ) ; 

E  =  getappdata (gcbf , ’E’ ) ; 
b  =  getappdata (gcbf , ‘b ') ; 
c  =  getappdata (gcbf ,  Vc ’) ; 

SF  =  getappdata (gcbf SF ’) ; 
eps  =  getappdata (gcbf , ’eps ') ; 

DkKnee  =  str2double (get (handles . editl4 , 'String' ) ) *1000000; 
dadnKnee  =  str2double (get (handles . editlS , ' String ' ) ) ; 

RRat  =  str2double (get (handles . edit 13 ,' String ')) ; 
kmx  =  DkKnee,/ (1 -RRat)  ;  %-kres; 
kmn  =  (DkKnee -kmx) /  %-kres; 

xnum  =  2 ; 
clear  RRat; 

[maxstrs,minstrs,minstrn,maxstrh]  =  rice (xnum, Ds tar , kmx, kmn, n, E, sy) ; 
maxstrs  =  maxstrs (1) ; 
minstrs  =  minstrs(l); 
maxstrn  =  maxstrn(l); 
minstrn  =  minstrn(l); 
dstrs  =  (maxstrs-minstrs) / 
dstrn  =  (maXstrn-minstrn) ;  . 

rb  =  minstrs/maxstrs ; 
prd  =  {dstrs*dstrn) / (2* (l-rb) ) ; 
clear  DkKnee; 
clear  minstrs; 
clear  maxstrs; 
clear  minstrs; 
clear  maxstrn; 
clear  minstrn; 
clear  dstrs; 
clear  dstrn; 
clear  rb; 
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%  TABLE  IS  FOR  Nf  from  le  -5  to  le  13,  (le'^  to  x)  ’’plastic  part",  (x  to 

le^'^)  "elastic  part". 

SM(1,2)  =  lel5; 

SMel  =  ({SF^2)/E)*((2*SM(l,2))^(2*b)) ; 

SMpl  =  (SF*eps)*{(2*SM{l,2))"(b+c)); 

SM(1,1)  =  SMel;  %  +  SMpl;  %  ELASTIC  PART 
SM(2,1)  =  prd; 

SM(2,2)  =  (Dstar/dadnKnee) ; 

SM(3,2)  =  le-5; 
clear  SMel; 
clear  SMpl; 

SMel  =  (  (SF^2) /E) * ( (2*SM{3,2) ) ^ (2*b) )  ; 

SMpl  =  (SF*eps) * ( (2*SM(3,2) ) ^ (b+c) ) ; 

SM(3,1)  =  SMpl;  %  +  SMel;  %  PLASTIC  PART 
clear  dadnKnee; 
setappdata(gcbE,'SM',SM); 
set (handles .edit29, ' string’ , SM {;, 1) ) ; 
set (handles.edit30, ’string' ,SM( : ,2) ) ; 
hold  on; 

axes (handles . axes3 ) ; 

%  Plotting  equation  and  the  three  point  solution  on  the  GUI- 
%  SMI  is  the  variable  with  100  point  calculated  equation 
%  SM  is  the  three  point  SWT  table  with  a  knee  point 

loglog(SM( : ,2) ,SM( : ,1) , ' :ok’ ) ; 

N  =  logspace (15, -5, 100)  '  ; 
lim  =  length (N) ; 
for  i  =  l:lim 

RHS  (i,  1)  =  (SF^2)  /E)  *  (  (2*N(i,  1)  )  ^  (2*b)  )  +  (SF*eps)  *.(  (2*N(i,  1)  )  ^  (b+c)  )  ; 
end 

SMl(l:lim,l)=RHS(l:lim,l); 

SMI (l:lim,2)=N(l:lim,l) ; 
loglog(SMl  (  :  ,2)  ,SM1  (  :  ,1)  ,  ’  -b’  )  ; 

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
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Main  subroutine  for  the  overall  calculations  and  control  of  the  program 


%================================================================== 

%  -  Executes  on  button  press  in  Runcalc. 

function  Runcalc_Callback (hObj ect ,  eventdata,  handles) 

%  hObject  handle  to  Runcalc  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 


w  =  str2double (get (handles .edits String ')) ;  %  Units  =  m 

af  =  str2double (get (handles . editl9, ' String' )) / 

Dstar  =  str2double  (get  (handles  .editl6,  'String' ))  ;  %  Units  ==  m 

setappdata (gcbf ,  'Dstar',  Dstar); 
sy  =  getappdata (gcbf ,  'sy'); 

n  =  getappdata (gcbf ,  'n'); 

ai  =  str2double (get (handles . editlO String ')) ; 

no  =  round ( (af-ai) /Dstar) ;  %  total  no.  of  elements  in  the 

specimen  rounded  off  to  nearest 
integer 

setappdata (gcbf ,  'no',  no); 

E  =  getappdata (gcbf ,  'E'); 

geom  =  get (handles . Li stbox_Geometry, 'value ') ; 

Idc  =  getappdata (gcbf ,' Idc ') ; 

Dkth  =  getappdata (gcbf,  'Dkth'); 
b  =  getappdata (gcbf , 'b' ) ; 
c  =  getappdata (gcbf ,' c ') ; 

tempxnum  =  str2double (get (handles . editSl ,' String ')) ;  %  from  GUI 

setappdata (gcbf ,' tempxnum' , tempxnum) ; 

DAM  =  zeros ( (no+300) , 1) ; 

DAMAGE  =  zeros(l,10); 

LNGTH  =  zeros (10,1); 

CYC  =  zeros (10, 1) ; 

KR  =  zeros ( (no+100) , 1) ;  %  KR  will  be  used  as  K  res 

setappdata (gcbf KR KR) ; 
retard  =  0; 
cnt  =  0; 

no_el_of f set  =  O; 
init  =  1; 

for  iterl  =  1 :no 

LNGTH (iterl, 1)  =  ai+  iterl*Dstar;  %  Storing  end  points  of  all 

elements  in  memory 

end 

ov  =  0; 

% - - - 

%  MAIN  FOR  LOOP  OF  CALCULATION  STARTS 

% - - - - - ^ - - 

for  p  =  l:no 

%  ovpercent  =  str2double (get (handles . editl7 ,' String ')) /lOO / 

%  determines  at  what  %  of  (af  ai)  to  apply  overload 
%  chk  =  round (ovpercent *no) ; 
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ovbrokenels 
ovcl  =  0; 


0/ 


THIS  IF  LOOP  APPLIES  OVERLOAD  (calls  the  Overload  subroutine) 


(V  ==  0  %  if  ov  is  1  then  overload  is  not  applied 

a  =  ai+(  (p-i) *Dstar) ; 
al  =  a/w; 
if  a  >=  0 . 0254 
clear  DeltaKt; 


[DAM, CYC, ovbrokenels, savedstrs,DeltaKt]  = 

°'^®^loa«3dam  ( a ,  al ,  p ,  DAM ,  no ,  LNGTH ,  Kmax ,  Kmin ,  CYC ) 

%  STORES  DAMAGE  FOR  OVERLOAD  1  CYCLE 


if  ovbrokenels  >=  1 

da  =  (ovbrokenels*Dstar) ; 
number  =  ovbrokenels-1; 
for  it4  =  p:p+number 

DADNA{it4,l)  =  a+da; 
DADMDK(it4,l)  =  DeltaKt/lOOOOOO; 
DADNDK(it4,2)  =  (da/1) / 
DADNA(it4,2)  =  DADNDK ( it4 , 2 ) / 

end 

end 


p  =  p+ovbrokenels; 
a  =  ai+ ( (p-i) *Dstar) ; 

residual  method  ############################ 

^  ffffff#ffff####Calls  the  weight  function  subroutine  ################# 

[K_res]  =  calculateKres (a,savedstrs,n,E,Dstar,sy) ; 

clear  Kres/ 
clear  XX; 
pb  =  p-i; 

for  it9  =  1;  length  (K__res) 

KR(it9+pb,l)  =  K_res (it9, 1) ; 

end 

setappdata (gcbf , ' KR ' , KR) ; 

clear  K_res; 


ovcl  =  1; 

OV  =  1; 

retard  =1; 

controls  if  loop  for  applying  the  retardation  by  shifting  range 
%cnt  =1;  3a 

end 

end 

% - - 

a  =  ai+ ( (p-i) *Dstar)  ; 
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al  =  a/w; 


clear  K_max; 
cledr  K_tnin; 
clear  applDK; 

if  Idc  ==  0 

K_max  =  str2double (get (handles .editl2, 'String' )) *1000000; 
K_min  =  str2double (get (handles . editll String ')) *1000000 ; 
setappdata  (gcbf ,  '  K__max '  ,  K_max)  ; 

setappdata (gcbf ,  ' K_min ' ,  K_min) ; 

elseif  Idc  ==1 

S_max  =  str2double (get (handles .editl2, 'String' )) *1000000; 
S_min  =  str2double (get (handles .editll, 'String' )) *1000000; 
DS  =  S_max-S_min; 

[DeltaK,DeltaKt]  =  getDeltaK (a, al , DS , Dkth, geom) ; 

RR  =  (S_min/S_max) ; 

K__tnax  =  (DeltaK/ (1-RR)  )  ; 

K_min  =  K^max  -  DeltaK; 
setappdata  (gcbf ,  'K__max',  K_max)  ; 

setappdata  (gcbf ,  '  K_tnin '  ,  K_min)  ; 

end 

applDK  =  K__max-K_min; 

K_max  ~  K_max  +  KR{p,l); 

K_min  =  K_min  +  KR(p,l); 
if  K_min  <0 
K_min  =0; 

end 

setappdata  (gcbf ,  ’  K_max '  ,  K__max)  ; 
setappdata (gcbf , ' K_min ' , K_min) ; 

rp  =  (1/ (  (1+n)  *pi)  )  *  (  (K_max/sy)  "^2)  ; 
clear  xnum; 

xnum  =  round (rp/Dstar) ; 

%  No.  of  elements  ahead  of  tip  within  Plastic  zone 
xO  =  Dstar; 


if  p  ==  230 

Stopppp  =1; 

end 

clear  Kmax; 
clear  Kmin; 

Kmax  =  K_max; 

Kmin  =  K__min; 

if  Kmin  <  0 
Kmin  =  0; 

end 

setappdata (gcbf , 'Kmax' ,Kmax) ; 
setappdata (gcbf , ' Kmin ' , Kmin) ; 
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Kmax  =  getappdata (gcbf , ' Kmax* ) / 
Kmin  =  getappdata (gcbf , ' Kmin ' ) ; 
R  =  Kmin/Kmax; 


DADNDK(p,l)  =  applDK/1000000; 

%  Storing  unmanipulated  DK  but  working  on  the  manip.  DK 


rp  =  2* (1/ ( (l+n) *pi) ) * ( (Kmax/sy) ^2) / 

%  CRACK  IS  ARRESTED  IF  THE  rp  <  Dstar 
if  rp  <  Dstar 

errordlg ( 'Crack  cannot  grow  further', 'No  crack  growth on ') / 
uiwait; 

end 

clear  xnum; 

xnum  =  round (rp/Dstar) ; 

%  No.  of  elements  ahead  of  tip  within  Plastic  zone 

setappdata (gcbf ,  'xnum',  xnum) ; 

xmin  =  Dstar; 

xmax  =  xnum*Dstar; 

xnu  =  xnum-l; 

xO  =  Dstar; 

EL  -  zeros (xnum, 2) ; 

% - - - 

na  =  n; 

% - - - - 

%  getting  the  maxstrs  minstrs  maxstrn  minstrn  from  RICE  +  Elastic 

% - ■ _ 

clear  rpmel; 
clear  rpcel; 

[maxstrsr,minstrsr,minstrnr,maxstrnr, Ingth, rpmel, rpcel]  = 

rice (xnum, xO, Kmax, Kmin, n, E, sy) ; 

% - - - . _ , _ 

xnum  =  length (maxstrsr) ; 
for  ita  =  l:xnum 

maxstrs (ita, 1)  =  maxstrsr (ita) ; 
minstrs (ita, 1)  =  minstrsr (ita) ; 
minstrn(ita,l)  =  minstrnr (ita) ; 
maxstrn(ita, 1)  =  maxstrnr (ita) ; 

end 

clear  maxstrsr; 
clear  minstrsr; 
clear  maxstrnr; 
clear  minstrnr; 

for  V  =  l:xnum 

EL(v,l)  =  maxstrs (v, 1) ; 

EL(v,2)  =  minstrs (v, 1) ; 

end 

SM  =  getappdata (gcbf ,  'SM')/ 

SF  =  getappdata (gcbf ,  'SF'); 

%  if  EL(1,1)  <  SF 


%  max  strs  associated  to  each  element 
%  min  strs  associated  to  each  element 
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for  it2  =l:xnum 

destrs  =  EL (it2 , 1) -EL (it2 , 2)  ; 

destrn  =  maxstrn (it2 , 1) -minstrn (it2 , 1) ; 

rbar  =  (EL (it2 , 2) /EL (it2 , 1)  )  ;  , 

■fit(it2,l)  =  (  (destrs*destrn)  /  (2*  (1-rbar)  ) )  ; 

%  Asigma  *  Aepsilon  prod  to  find  Nf 

end 

clear  it2; 

life  =  gettingNf (f it) ; 

%  Interpolated  values  of  Life  for  the  'fit'.. on  SM  table  curve 

ni  =  (l-DM  (p,  1)  )  *life  (1, 1)  ;  %  ni  for  current  element  to  fail 

pa  -  p-l;  %  just  subtracting  1  from  p 

if  ni  <  1 
ni  =  1; 

end 

It  =  length (life) ; 
for  it2  =l:lt 

DAM(  (pa+it2)  ,  1)  =  (ni/lif  e  (it2 , 1 )  )  +  DAM  (  (pa-f  it2 )  ,  1)  ; 

%  stores  the  damage  value  of  element  upto  current  element  broken 
end 

clear  it2; 

brokenels  =  length (find (DAM (p :  (pa+lt) )  >=  0.99999999)); 

%  Finding  the  no.  of  broken  elements  to  calculate  the  crak  growth  rate 

%  for  it3=p:no  ' 

%  ■  CYC (its, 1)  =  ni  +  CYC(it3,l); 

%  Stores  the  total  no  of  cycles  for  element  p  at  or  after  breaking 
%  end 

number  =  brokenels -1; 
if  number  <0 
number  =  0 ; 

end  ‘ 

da  =  (Dstar  +  (number) *Dstar) ;  . 

.  DelK  =  DADNDK(p, 1) ; 
for  it  =  p:p+number 

DADNA(it,l)  :=  a+da; 

DADNDK(it,l)  =  DelK; 

DADNDK(it,2)  =  (da/ (ni+ovcl) ) ; 

DADNA(it,2)  =  DADNDK (it , 2 ) ; 

end 

p  =  p+number; 
ovc 1=0; 

P  =  P+1; 
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if  p  ==  round (init*no/100) ; 

string  =  ['calculation  is  ' , num2str (init ) , ' %  complete']; 
wb  =  waitbar {0 , string)  ; 
waitbar (p/no) ; 
close (wb) ; 
init  =  init+1; 
if  init  ==  95/ 
hhh=l; 

end 

end 


end 

setappdata (gcbf , 'DADNA' ,  DADNA) ; 
setappdata (gcbf , ’ DADNDK ' ,  DADNDK) ; 
hold  on; 

DADNDK  =  get appdata (gcbf , 'DADNDK' ) ; 
axes (handles . axesl) ; 


if  Idc  ==  1 

loglog (DADNDK ( (1 : round (no/20) :no) ,1) , DADNDK ( (1 ; round (no/20) :no) ,2) , ' - 

.*r’) 

set (handles . edit20 , ' string ' , DADNDK (1 :no, 1) ) ; 
set (handles . edit21, ' string ' , DADNDK (1 :no, 2) ) ; 

elseif  Idc  ==  0 

plot (DADNA ( (1 : round (no/20)  :no) ,1)*1000, DADNA ( (1 ; round  (no/20)  :no)  ,2)  ,  '  - 

*b'  ) 

set (handles. edit20, 'string' ,DADNA(l:no, 1) *1000) ; 
set (handles .edit21, ' string' ,DADNA(l:no, 2) ) ; 

end 

setappdata (gcbf ,' count ' ,  1)  ; 


h+++++++++++++++++++++H 


+++++-f+++ +++  +++++++++  + 


Subroutine  for  the  application  of  the  overload  or  block  load  change 


%  THIS  FUNCTION  CALCULATES  THE  DAMAGE  DUE  TO  OVERLOAD 


function [DAM, CYC, ovbrokenels,savedstrs,DeltaKt]  = 

overloaddam (a, al,p,DAM, no, LNGTH, Kmax, Kmin, CYC) ; 

E  =  getappdata (gcbf ,  'E*)/ 

sy  =  getappdata (gcbf ,  'sy')/ 

n  =  getappdata (gcbf ,  ’n’)/  . 

SF  =  getappdata (gcbf ,  ’SF')/ 

eps  =  getappdata (gcbf ,  'eps'); 

Dstar  =  getappdata (gcbf , 'Dstar' ) ; 

Dkth  =  getappdata (gcbf , 'Dkth' ) ; 

SM  =  getappdata (gcbf ,  'SM'); 

ofdkmax  =  Kmax; 
oldkmin  =  Kmin; 

Kmax  =  Kmax*2; 

Kmin  =  Kmin; 

R  =  Kmin/Kmax; 

ep  =  getappdata (gcbf ,  'ep');  ; 

DeltaKt  =  Kmax-Kmin; 

%  >  - - - 

rp  =  2*(l/((l+n)*pi))*((Kmax/sy)^2); 

if  rp  <  Dstar 

errordlg (’ Crack  cannot  grow  further', 'No  crack  growth on ') ; 
uiwait; 

end 

clear  xnum; 

xnum  =  round (rp/Ds tar) ; 

%  No.  of  elements  ahead  of  tip  within  Plastic  zone 

setappdata (gcbf ,  'xnum',  xnum) ; 

xmin  =  Dstar; 

xmax  =  xnum*Dstar; 

xnu  =  xnum-1; 

xO  =  Dstar; 

EL  =  zeros (xnum, 2 ) ; 

b  =  getappdata (gcbf , 'b' ) ; 
c  =  getappdata (gcbf ,' c ') ; 
na  =  n; 

% - - - - - - - - 

%  getting  the  maxstrs  minstrs  maxstrn  minstrn  from  RICE  +  Elastic 
- - - - - - - 

[maxstrsr , minstrsr , minstrnr , maxstrnr]  =  rice (xnum, xO, Kmax, Kmin,n, E, sy) ; 

for  ita  =  l:xnum 

maxstrs (ita, 1)  =  maxstrsr (ita) ; 
minstrs (ita, 1)  =  minstrsr (ita) ; 
minstrn (ita, 1)  =  minstrnr (ita) ; 
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maxstrn (ita, 1)  =  maxstrnr (ita) ; 

end 

clear  maxstrsr; 
clear  minstrsr; 
clear  maxstrnr; 
clear  minstrnr; 

savedstrs  =  zeros (xnum, 1) ; 
savedstrs ( : , 1)  =  minstrs ( : , 1) ; 

EL(l:xnum,l)  =  maxstrs (1 :xnum, 1) ; 

%  max  strs  associated  to  each  element 
EL(l:xnum,2)  =  minstrs (1 :xnum, 1) ; 

%  min  strs  associated  to  each  element 

for  it2  =l:xnum 

destrs  =  EL (it2 , 1) -EL (it2 , 2)  ; 

destrn  =  maxstrn (it2 , 1) -minstrn (it2 , 1) ; 

rbar  =  (EL {it2 , 2 ) /EL (it2 , 1) ) ;  ‘ 

fit(it2,l)  =  ( (destrs*destrn) / (2* (1-rbar) ) ) ; 

%  Asigma  *  Aepsilon  prod  to  find  Nf 

end 

life  =  gettingNf (f it ) ;  %  Getting  Nf  values  from  the  SM  table 

ni  =  1;  %  ni  =  1  for  one  cycle  overload 

pa  =  p-l;  %  just  subtracting  1  from  p 
for  it2  =p: (xnu+p) 

DAM(it2,l)  =  (ni/life (it2-pa, 1) ) +  DAM{it2, 1) ; 

%  stores  the  damage  value  of  element  upto  current  broken 
end 

%  CHECKS  FOR  HOW  MANY  ELEMENTS  THE  DAMAGE  GOES  ABOVE  1 
ovbrokeneis  =  length (find (DAM (p : (xnu+p) ) >=0 . 99999999) ) ; 

%  for  itB=p:no 

%  CYC(it3,l)  =  ni  +  CYC{it3,l);  %  Stores  the  total  no  of  cycles 

for  element  p  at  or  after  breaking 
%  end 

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
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Subroutine  to  calculate  stress  profiles  ahead  of  the  crack  using  Rice  and 
Elastic  solution. 


%  RICE  AND  ELASTIC  SOLUTION  SUBROUTINE 


function  [maxstrs , minstrs, minstrn, maxstrn, Ingth, rpmel, rpcel]  = 

rice (xnum,xO; Kmax,Kmin,h,E, sy) ; 


ey  =  sy/E; 


if  xnum  <=  1 
xnum  =10; 

end  ’  • 

for  i  =  l:xnum 

lngth(i,l)  =  i*xO; 

end 

pow  =  n/ (1+n) ; 
spow  =  1/ (1+n) ; 

for  it  =  l:xnum 

S(it,l)  =  sy* ( (Kmax*2) / ( (1+n) *pi* (sy*2) *lngth (it, 1) ) ) ^pow; 
EP(it,l)  =  ey* ( (Kmax^2) / ( (1+n) *pi* (sy^2) *lngth(it, 1) ) ) ^spow; 

end 

rpmel  =  length (find (S (: ,1)  >=  sy) ) ; 
if  rpmel  ==  0 

xm  =  ( (Kmax/sy) ^2) * (1/pi) * ( (1/ (1+n) ) -0 . 5) ; 

%offset  for  the  elastic  solution  to  be  continuous  at  r^  . 
else 

for  itl  =  1: rpmel 

maxstrs (itl)  =  S(itl,l); 
maxstrn (itl)  =  EP(itl,l); 
end 

xm  =  ( (Kmax/sy) ^2) * (1/pi) *( (1/ (1+n) ) -0 . 5) ; 

end 


%  elastic  solution  passing  through  same  point  has  to  be  found 
for  it2  =  rpmel+l:xnum 

maxstrs (it2)  =  Kmax/ (sqrt (2*pi* (Ingth ( it2 , 1) - (xm) ))) ; 

maxstrn (it2)  =  maxstrs ( it2) /E; 

end 


%  FINDING  DELTA  SOLUTION 

% - - - - - - 

DK  =  Kmax-Kmin; 

for  it3  =  l:xnum 

DS  (it3, 1)  =  2*sy* ( (DK^2) / (4* (1+n) *pi* (sy^2) *lngth (it3 , 1) ) ) ^pow; 

DEP  (it3 ,1)  =  2*ey* ( (DK^2) / (4* (1+n) *pi* (sy^2) *lngth (it3 , 1) ) ) ^spow; 

end 
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rpcel  =  length (find (DS (:, 1)  >=  2*sy) ) ; 
if  rpcel  ==  0 

xc  =  ((DK/sy)^2)*(l/(4*pi))*{(l/(l+n))-0.5) ; 

%offset  for  the  elastic  solution  to  be  continuous  at  rc 
else 

for  it4  =  1: rpcel 

deltastrs  {it4)  =  DS(it4,l); 
deltastrn  (it4)  =  DEP(it4,l); 
end 

xc  =  (  (DK/sy)  *^2)  *  (1/ (4*pi)  )  *  (  (1/ (1+n)  ) -0  .5)  ; 

end 

%  elastic  solution  passing  through  same  point  has  to  be  found 
for  its  =  rpcel+l:xnum 

deltastrs  (its)  =  DK/ (sqrt (2*pi* (Ingth (its , 1) - (xc) ) ) ) / 
deltastrn  (its)  =  deltastrs (its) /E; 

end 

% - - - - - - - 

%  MINIMUM  STRESS  PROFILE  and  MINIMUM  STRAIN  PROFILE 

% - - - - - 

for  its  =  Irxnum 

, minstrs (its)  =  (maxstrs (its)  -  deltastrs  (its) )  ;  ‘ 
minstrn(itS)  =  (maxstrn (itS)  -  deltastrn (its) ) ; 

end 
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Subroutine  to  calculate  the  Residual  SIF  associated  with  the  elements 
using  the  Weight  Function  solution 


%  This  function  takes  input  from  the  Overload  affecting  loop  to 
%  calculate  the  residual  stress  difference  profile,  and  then  goes 
%  on  to  calculate  the  weight  function  associated  with  it 


function  [K___res]  ==  calculateKres  (a,  savedstrs ,  n, E, Dstar,  sy)  / 
numm  =  length (savedstrs) / 
for  it6  =  IrnumTTt; 

LL(it6,l)  =  a  +  (Dstar*it6);  %  end  points  of  the  elements 
end 

clear  it6/ 
crk_ln  =:  a; 

YY{1,1)  =  savedstrs (1, 1) ; 

XX(1,1)  =  a; 
for  it  =  2:numm+l 

YY(it,l)  =  savedstrs ( (it-1) , 1)  ; 

XX(it,l)  =  LL(it-l)  / 

end 

Int  =  length  (XX)/ 
for  it2  =  2 : Int 

%  it2  represents  the  position  of  the  imaginary  crack  tip 
clear  u; 
clear  crk; 
clear  Ml; 
clear  M2; 
clear  M3; 


crk  =  XX {it2 , 1) ; 

for  its  =  l:it2 

u(it3)=XX(it3,l)/crk; 

end 

clear  it3; 

K_temp  =  0; 
alpha  =  (crk/0 .5) ; 
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%  Calculation  of  Ml, M2  and  M3  for  the  edge  crack  plate 

Mi  =  ((-0.029207+alpha*(0.213074+alpha*(-3.029553+alpha*(5. 901933+ 
alpha* (-2.657820) ))))  / 

(1+alpha* (-1.259723+alpha* (-0 . 048475+alpha* ( . 48125+alpha*  (- 
.526796+alpha* ( .345012)  ))))))  ; 

M2  =  (  ( . 451116+alpha* (3 . 462425+alpha* (-1 . 078459+alpha* (3 . 558573  + 

alpha* (-7.553533) ))) )/ 

(1+alpha* (-1 . 496612+alpha* ( . 764586+alpha* (- . 659316+alpha*  ( 
0.258506+alpha*(. 114568))))))  ); 

M3  =  (  (.427195+alpha*(-3.730114+alpha*(16.276333+alpha*(- 
18 . 799956+alpha* (14 . 112118) ))))  / 

(1+alpha* (-1 . 129189+alpha* ( . 033758+alpha* ( . 192114+alpha*  (- 
0.658242+alpha* (.554666) )))) )  ); 


FI  =  0; 

■  F2  =  0; 

for  in=l : (it2-l) 

A=  ((YY(in+l,l)-Yy(in,l))/(u(in+l)-u(in))) ; 

B=  yY(in,l)  -  (A*u(in)  )  ; 
test  =  A*crk+B; 

Fl=  (crk/  ( (2*pi*crk)  .^0.5)  )  .*  (A*  (  (-4*  (l-u(in+l) )  .''0.5)  -  (2*M1*  (1- 
u(in+l) ) )  +  (4/3* (l-u{in+l) )  .*1.5) - (4/3*M2* (1-u (in+1) )  . *1.5)  .  .  . 

+ (Ml* (1-u (in+1) ) . *2) - (M3* {l-u(in+l) ) . *2) + (4/5*M2* (1- 
u(in+l) ) .*2.5) + (2/3*M3* (1-u (in+1) ) . *3) ) -B* ( (4* (1-u (in+1) ) .*0.5) + 
(2*M1* (1-u (in+1) ) )  +  (4/3*M2* (1-u (in+1) )  . *1 . 5)  +  (M3* (1-u (in+1) )  .  *2) ) ) 

F2=  (crk/ ( (2*pi*crk) .*0.5) ) .* (A* ( (-4* (l-u(in) ) .*0.5)- (2*M1* (1- 
u (in) ) ) + (4/3* (l-u(in) ) .*1.5) - (4/3*M2* (1-u (in) ) . *1.5) . . . 

+  (Ml*  (1-u  (in) )  .  *2)  -  (M3*  (1-u  (in) )  .  *2)  +  (4/5*M2*  (1- 

u(in)) .*2.5)+(2/3*M3* (l-u(in) ) . *3) ) -B* ( (4* (1-u (in) ) .*0.5) .. . 

+(2*M1* (1-u (in) ) )+(4/3*M2*(l-u(in) ) . *1 . 5 ) + (M3* (1-u (in) ) .*2) ) ) ; 

if  ((isnan(Fl))  |  (isnan(F2))) 

FI  =  0; 

F2  =  0; 

end 

K___temp  (in)  =F1-F2  ; 

end 

%  total  stress  intensity  factor  is 
alt  =  it2-l; 

K_res(alt,l)  =  sum(K_temp); 
clear  K_temp; 
clear  alt; 

end 
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Subroutine  to  calculate  the  SIF  solution  for  the  applied  loading 


%  RETURNS  THE  DeltaK  VALUES  BACK  TO  THE  CALL 
function [DeltaK,DeltaKt]  =  getDeltaK(a,al,DS,Dkth,geom) 
switch  geom 

case  1  %  CENTRE  CRACK 

F=  (l-(0.5*al)+(0.326*(al^2)))/sqrt(l-al) ; 

%  F  Solution  from  Norman  Dowling 

DeltaKt  =  F*DS*sqrt {pi*a) ; 

case  2  %  EDGE  CRACK 

F=  (0.265*(l-al)^4)+((0.857+0.265*al)/(l-al)^1.5) ; 

%  F  solution  from  Norman  Dowling 
DeltaKt  =  F*DS*sqrt (pi*a) ; 

case  3  %  DOUBLE  EDGE  CRACK 

F  =  1+ (0.122* (cos (pi*al/4) )^4) )* (sqrt( (2/(pi*al) ) *tan  (pi*al/2 ) ) ) ; 
%  F  solution  from  Norman  Dowling 

DeltaKt  =  F*DS*sqrt (pi*a)  ; 

end 

%.: _ _ _ _ _ _ 

if  DeltaKt  <=  Dkth 

err  =  'Delta  K  is  less  than  threshold  value  =  '; 

Delt  =  num2str (DeltaKt/1000000, 3)  ; 
errd  =  strcat (err, Delt) ; 
errordlg(errd,'Applicationhalted','on'); 
uiwait; 

end 

%  If  any  adjustment  or  correction  is  needed  it  can  be  applied  to  this 
%  extra  variable  "DeltaK"  here  which  is  equal  to  DeltaKt  normally 
DeltaK  =  DeltaKt; 


